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SIMULATION OF THE FOKKER-PLANCK EQUATION BY RANDOM WALKS 
OF TEST PARTICLES IN VELOCITY SPACE WITH APPLICATION 


TO MAGNETIC MIRROR SYSTEMS 
by Gerald W. Englert 
Lewis Research Center 

SUMMARY 

The structure of an analytical expression obtained from a random-walk process was 
compared with that of the Fokker-Planck equation. The step sizes and probabilities of 
taking steps in various directions were related to the coefficients of the Fokker-Planck 
equation. Spherical coordinates with azimuthal symmetry about the polar axis were 
used. 

The technique was applied to a magnetic-mirror system for confining charged par- 
ticles. The cases selected have approximate analytical solution via the Fokker-Planck 
equation. Velocity distributions inside and end losses from this system were determined 
by the random-walk procedure and compared with results of the analytical solutions. 

The effect of reducing the maximum impact parameter cutoff distance well below the 
Debye length was studied. This reduction increased step sizes and decreased collision 
frequency in such a manner that the probability of a test particle walking to specified 
locations in velocity space remained approximately constant. A sampling of random 
numbers from the lower impact-parameter range greatly reduced computing time and 
yet permitted reliable distributions within 10- to 20-percent accuracy. 

The random-walk results were, in turn, permitted an appraisal of some simplifying 
assumptions made in the analytical solution of the Fokker-Planck equation. An iteration 
procedure was used to find a self-consistent solution of the distribution of test particles 
walking through an ensemble of field particles, on which the Fokker-Planck coefficients 
were based. Analytical solutions assume that the distribution of test particles and loss 
rates are insensitive to assumed field-particle distributions. For the cases studied this 
was found to be a good assumption. Two steps in the iteration process made the 
random-walk solutions self-consistent within the scattering of the data for a sampling 
size of 500 test particles. The popular assumption of separability of the distribution 
function was also verified. 

The random-walk method should apply with little additional difficulty to a number of 
problems that cannot be solved analytically. 



INTRODUCTION 


The Fokker- Planck equation, in general, describes the time development of a 
Markov process. Such a process is characteristic of the nature of classical collisions 
where each event depends on the present conditions and is independent of the past 
(ref. 1, p. 157, and ref. 2, p. 369). Written in velocity space, this equation provides 
a popular approximation to the collisional terms in the Boltzmann equation (refs. 3 to 5). 
Included in this approximation are first and second moments of velocity increments de- 
scribing two-body encounters (p. 75 of ref. 6). These moments enter as coefficients 
of the Fokker-Planck equation and permit it to describe changes in the velocity distri- 
bution function resulting from influence of dynamic friction and dispersion (refs. 

7 and 8). In general, the velocity distribution function is used to weight the moments 
of velocity increments, making the Fokker-Planck equation nonlinear and very difficult 
to solve for many problems of interest. 

Another approach to collision problems is through the study of random walks of 
test particles. This approach replaces the mathematical complexity with a relatively 
simple but repetitious process, which is greatly facilitated by high-speed computers. 
Usually the walks of a large sample of test particles through an ensemble of field par- 
ticles must be studied to determine how they distribute in velocity space. For a self- 
consistent solution, the resulting test-particle distribution should coincide with the 
field-particle distribution, requiring an iteration process. 

A random walk will be defined herein as a process in which step sizes and proba- 
bilities of taking various steps are average values (depending on field-particle distribu- 
tions and test-particle location in velocity space for the case at hand). In contrast to 
this procedure is the more popular Monte Carlo method. This technique makes random 
selections from distributions of pertinent variables (such as impact parameter and rela- 
tive velocity) at each collision and performs an interaction calculation at each step to 
determine the resulting test-particle location (refs. 9 and 10). Such a detailed calcula- 
tion at each step is excellent for physical and mathematical clarity, but would involve 
an excessive amount of computer time for cases involving a very large number of en- 
counters. On the other hand some physical clarity can be easily lost in the representa- 
tion by random walks, especially in determining the probability of taking steps in vari- 
ous directions. The ability to relate step sizes and probabilities to the well explored 
coefficients of the Fokker-Planck equation would thus be very useful. 

The random-walk and the Fokker-Planck concepts depend primarily on the same 
combinatory laws of probability; however, the random walk as depicted herein is re- 
stricted to steps on a grid. A detailed generation of the Fokker-Planck equation from a 
random-walk procedure is available in the literature only for the one-dimensional prob- 
lem with constant coefficients. Reference 2 (ch. 14) shows that for this case the proper 
step size of a random walk is dependent on the second moment and that the probabilities 
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of taking steps in positive and negative directions depend on both first and second mo- 
ments. This derivation is extended herein to problems involving three dimensions and 
variable coefficients. In this more general case the step sizes and probabilities are 
again determined in terms of Fokker-Planck coefficients. The Fokker-Planck coeffi- 
cients are evaluated for the binary Coulomb encounters used in plasma physics. 

A main difficulty in obtaining final results is due to the long range nature of Coulomb 

encounters, which makes the step sizes of the walks extremely small. For example, at 

7 

typical experimental or laboratory conditions, a particle may average 10 encounters 

(steps) before it reaches a 90° deflection from its initial direction. At thermonuclear 

1 fi 

conditions this number easily reaches 10 steps. Tracing such a large number of 
steps for a representative number of test particles would require an excessive amount 
of time even for today's fastest computers. A selective sampling technique of some 
nature must be therefore employed to the random-walk procedure. 

The coefficients of the Fokker-Planck equation are insensitive to the cutoff assumed 
for the collisional impact parameter (refs. 11 and 12). As the maximum impact param- 
eter is reduced, the average step size increases, but the number of encounters (steps) 
per unit time decreases. With a low impact-parameter cutoff distance a test particle 
would take a relatively small number of large steps. On the average, however, it would 
reach a specified boundary in velocity space in about the same time as in a walk based 
on a higher impact-parameter cutoff and requiring a large number of small steps. The 
computing time, however, varies inversely as about the square of the step size and is 
thus much less for walks based on lower impact-parameter cutoff distances. 

To study results of this sampling technique, it was applied to the end-loss problem 
of magnetic -mirror systems. The scattering of charged particles into a loss cone has 
received considerable attention over a period of years. However, because of the com- 
plexity of the Fokker-Planck equation, the rate of such loss has been determined ana- 
lytically for only the simplest of cases (refs. 5 and 3). Numerical solutions (ref. 13) 
are available for only a restrictive set of initial and boundary conditions. Desired also 
are particle distributions inside the mirror system for use in stability studies (ref. 14). 

In summarizing, the attempt herein is to exchange a difficult analytical problem 
with a fairly straightforward computing procedure. Computing time is made reasonable 
by a sampling technique. Considerable physical detail can be incorporated into such a 
procedure with little increase of mathematical complexity. 


ANALYSIS 

The general form of a random-walk equation will be derived for walking on a co- 
ordinate grid in N dimensions. (One-dimensional walks with constant step sizes and 
probabilities are treated in refs. 2 and 8. Ref. 2 writes the results in the form of a 



Fokker- Planck differential equation. Ref. 8 uses the random walk to derive solutions 
to some boundary value problems, but does not arrive at a differential equation as an 
intermediate step. The Fokker-Planck equation is also derived in ref. 8, but not from 
the concept of a random walk on a grid. ) It is not necessary for the grid to be orthog- 
ornal or have constant spacing. Comparison will then be made with the Fokker-Planck 
equation applied to inverse square collisions in velocity space. The results will be 
adapted to a study of charged particles in a magnetic-mirror system. 

A list of symbols is given in appendix A. The International System of Units (SI) 
will be used throughout with the exception of temperature T being reported in keV and 
the corresponding Boltzmann constant k in joules per keV. 


Development of the Random-Walk Equation 

Let • • •» independent coordinates that span the space of interest. 

Let p t and q t be the probabilities that the i th coordinate will increase or decrease 
during the course of a step. To identify these conditional probabilities, only their loca- 
tion at the start of a step will be labeled in their arguments. For example, compared 
with the more conventional nomenclature, 


Pl(£j> ^2’ ' ■ ■ > 


means 


Pl^l + ' ’ ’’ ^2’ ' ‘ ' ’ 


and 


%2’ ' ' ' 5 


means 


^ 2 ’ ' ' * > ^ 2 ’ ' ' ■ ’ 
it. 

Let the probability of a step being in the ±i u direction be 1/N; that is, 

Pj(^l> ^2’ ■ ■ ■ ’ ^i^ 1’ ^2’ ’ ’ •» = ~ (1) 

Le ‘ V«l> «2 «N> denote the chance of a particle being at location 
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(a) Steps along component direc- 
tions; number of components 
changed per encounter, 1. 



Figure 1. - Examples of random walks in two- 
dimensional Cartesian coordinates. 


(Ij, £ 2 , . . . , £ n ) at the start of the n th step and At be the time interval between 
starts of the n and n + v steps. 

Consider the case for tj = 1, illustrated in figure 1(a). The change in the number of 
particles at (Ij, $ 2 > • • • > l N ) during At is the probability of reaching (|p l 2 , . . • , 
£ n ) from the closest neighbors minus the chance of a particle leaving (Ip l 2 , • • ■ » £jsf) 
and going to a nearest neighbor. Since a step of finite size is taken each time increment, 
there is no chance of a particle standing still. The conservation of particles at (ijj, 
l 2 , . . . , l N ) can thus be written 
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u m-l^l’ h’ ‘ 

• • ’ ^ “ u n^l’ h’ ‘ 

• • > 4^) ~ - A ^i> 4 2 > ' 

• • ’ ^ u n^l " A ^l’ h’ • • 



+ +^i> $2’ ' ' 

' * ^N^ u n^l + A ^r ^2’ ‘ ’ 

•» 4 N ) 



+ P2^1> ^2 ” A ^2’ ' ‘ 

• > ^N^n^l’ h ' A ^2> • ' 

•* 



+ ^2^1’ ^2 + A ^2’ ’ ‘ 

•» ^ u n^l’ h + A h’ ■ • 

•» %) 



+ . . . +P N Up t 2 > * * *’ ~ A ^ u n^V h’ * * *’ 


+ q N (^i» %2’ • ■ 

1 ‘ ’ ^N + A ^N^ u n^l> 

4 2 ’ • • •> 4 N + a 4 n ) 




- Pi(^i> %2’ ■ ■ 

■ > ^V^l’ ^2’ • • 

■ ’ i> 4 2 > • • 


4 2 . • • 

• > 4^ 

- P2^1> %2’ ■ ■ 

•’ ^N^n^l’ h’ ■ ■ 

‘ ’ 4j^) " ^2^1* %2’ ‘ ‘ 

•» 4^ u n :^1* 

4 2 » • • 

•> ^ 

- • • • - P N (4 1; 

l 2 > • • •> ^ u n^l’ 

h’ • ‘ • • ^ " q N^i’ 

4 2 > ••••<: 

N^ u n^i> 

z 2 > • • 


where use is made of the first and second laws of composition of probability. 
Multiplying through by 2N and adding and subtracting u n terms results in 

2N t u n+l^l 5 ^2’ * * •> In^ " u n^V ^2’ ' • ’ ’ In^ = Np l^l ” A lp l 2 ’ * * *’ ^ u n^l " A lp l 2 ’ * ’ ’ ’ 

- [1 - N Pl U x - A| r * 2 , . . ^ N )u n (^ - A* p | 2 | N ) 

+ Nq 1 (| 1 + A| 1 , | 2 , . . Ij^yii + A lp l 2 ’ * * * ’ 1^ 

- [1 - Nq x (^! + A| p | 2 , . . + A*i, l 2 > * • - *N> 

+ u n^l“ A ^P ‘ ‘ 2u n^l 5 l 2 ’ ’ * *’ W + u n^l + A lp l 2 ’ * * * ’ 1^ 

+ Np 2 fep - A| 2 , . . I^yij* l 2 “ A ^2’ ’ * ’ ’ W 

- [ 1 - Np 2 (£ 1 , l 2 " A y * ’ *’ yK^l’ ^2 ” A ^2’ ’ * *’ 1^ 

+ Nq 2 (4p l 2 + Al 2 > * • •> lN^ u n^r l 2 ’ * ‘ ' * W 

-[1 - Nq 2 (| 1 , ^ 2 + A^ 2 , ■ ^N , lJ u n^l> l 2 + A ^2" ' * * > 

+ u n (ll> l 2 " A ^2 7 * * ' * In^ 2 y^l» h’ ' ’ ' » + u n^l» ^2 + A ^2’ * ' * » 1^ 

+ . . , + Npjjdj, l 2 » * • • > In “ A %^ u n^l» l 2 » ■ ■ ’ » In ” a In^ 

- f 1 - n Pn(Ii» l 2 ’ * * * * In "* A lN)] u n(lr l 2 > - ■ ~ a In^ 

+ N yir l 2 > * • •’ «N + Wi> « 2 ’ • • •’ !n + a !n> 

- [! - Nq N (|j, l 2 , - * • , l N + A l N )]yip l 2 * * * •» In + A In^ 

+ yip l 2 » ‘ In " a In^ " 2 yii» l 2 » * * •» l N ) + yir l 2 » * • * j In + a In^ 
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Using equation (1), combining p^ and terms of like arguments and i subscripts, 
collecting u n terms and finally writing the result in finite difference notation yields 


2 A^ An = 2 M At= ‘ 2 { 6l[(Pl ’ qi)Un]+62[(P2 ’ q2)U2]+ ' ‘ • +6 N [(p N' C lN )u j} + ^ (6 ll U u^22 u n + - ' ’ + 6 NN U n } 
which approaches the limit 


N 

9u_ V j 9 [(Pj - qj> u ] A ^j t 1 a 2 u 
at 3«i A t 2 0|2 N At 

as step size (grid spacing) becomes arbitrarily small. 

This is the more conventional approach to a random-walk expression. However, 
next consider the case where steps are taken in groups of rj = N; one along each coor- 
dinate direction. (This corresponds to the later physical description of a test particle 
in polar coordinates where its magnitude and direction (v, 0, and cp) change each en- 
counter. ) This corresponds to steps across the diagonals of N-dimensional lattices 
which will herein be called path lines. This is illustrated for the case of N = 2 on 
figure 1(b). 

Repeating this procedure with 77 = N gives 

3u f'f At 1-la VAtp* 

at / A 3|. At 2 at 2 At 

iTf l 1 d k 

This same result could be obtained by letting 

Pj(£ ^2’ ‘ ' ‘ ’ + l’ ^2’ ‘ ’ ‘ 1 i = 1, 2, . . . , N 

in the derivation regardless of N, and letting r) = 1; that is, there is a sure chance of 
stepping either in the plus or minus direction of each of the components for each time 
increment At. 

Equation (2) is in the form of the Fokker- Planck equation except that no mixed 
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partial derivatives are present. Generation of these terms does not appear possible 
when the walks are restricted to the nearest neighbors on a grid. Their magnitude will 
be discussed later. 

•Hi 

Note that p. - is like a bias in the i in direction and that the terms involving 
first-order derivatives could be, for example, related to friction. To keep the process 
from being degenerate, p^ - q. must approach zero as step size approaches zero (p. 324 
of ref. 2). The terms involving second derivatives are similar to diffusion expressions 
with (A£.) /At similar to a diffusion coefficient. 

When a stochastic process is being simulated by a random walk on a grid, the step 
sizes are often statistical averages. It can be deduced from references 2, 8, 15, or 16 
that second moments determine the proper grid spacing. Consider, for example, the 
application discussed in the next section. In any small region of velocity space the test 
particle taking a random walk encounters many other field or background particles, one 
at a time. A field particle can be in any accessible region of velocity space as it starts 
its encounter with a test particle (its probable location is defined by appropriate distri- 
bution functions). The grid spacings are thus based on statistical averages of encounters 
with field particles (fig. 2), these averages nevertheless being dependent on test- 
particle location. 


pj and qj Probabilities of increase and decrease, 
respectively, of i th component 
(i = 1, 2, 3) 

^/(AVj) 2 component of step size 


v 3 
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Fokker-Planck Equation 


The widely accepted Fokker-Planck equation for studying plasmas in N dimen- 
sional Cartesian coordinates is 

N 

Z - — (f( Av.) ) +- — (f( Av. Av.) ) 

av. 1 2 av. av. 1 ] 

i,]=i J 

where (3f/at) c denotes the change of probability density (distribution function of test 
particles) with respect to time due to collisions. The coefficients (AVj) and (AVj Avj) 
are averages over distributions of field particles and scattering angle (ref. 11). The 
angular brackets signify a time average; that is, (Av-) is the average increment per 

± -L. 1 

unit time of the i Ln component of velocity. The average value per collision, denoted 
by a bar, is equal to the time average divided by the collision frequency v; for example, 
Av t = (AvJ/v. 

The Fokker-Planck equation is transformed to spherical coordinates with symmetry 
about the polar axis in appendix B. This is essentially the same coordinate system as 
used in parts IV and V of reference 11. In reference 11, cos 6 is used in place of the 
polar angle 9 used herein. The assumption of aximuthal symmetry is not too restric- 
tive for many problems in plasma physics and serves mainly to reduce cumbersome ex- 
pressions. The analysis follows in a straightforward manner without it. In this coor- 
dinate system 



/an 

\3tJ 


C V 


2 3v 


sin 0 30 


- A A (v 2 f( AV» - J_ A (f sin 0< A0» + -A + I — 9 2 f sin 9((Ag)j|) 

2v 2 av 2 2 sin 9 ae 2 

+ 1 a 2 v 2 f sin 9(Av A 8) _ JL_ _3_ |v 3 f[( (A6) 2 ) + sin 2 0( (a<^) 2 )]| 

+ v 2 sin 0 3v 39 2v 2 9v 

1 — / f sin 0[~2(A0 Av) - v sin 0 cos 0((A<p) 2 )]]} 

3in 0 30 \ J 


2v sin < 
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Using the identities 


3 2 v 2 f((Av) 2 ) = ((av) 2 / a 2 v 2 f + 2 _3_ / 2 f 3((Av) 2 )\ v 2 { e 2 <(Av) 2 ) 

av 2 av 2 / av 2 


and 


= <(Ae) 2 > a2 i sln 6 * 2 A (t sl „ e . t sl „ « a2 <^) 2 > 

90 2 90 2 a6> \ dd J dQ 2 


the Fokker-Planck equation can be written as 


m + i(. 

\St) c 2 V 


3 2 (( Av) 2 ) | 3 2 ((A 0 ) 2 ) 


3v 


dB‘ 


= - — — fv 2 f (Av) - ^ +- ((A 0 ) 2 ) +- sin 2 0 <(Aip) 2 )] 

/ y 2 3 v L Sv 2 2 J 


-^4s in *U) 

sin 0 30 L \ d9 


_ (A 0 Av) t sin 9 cos 0 


1 3 2 v 2 f Sin d( Av A 0 ) 


o 

v sin 0 


3 v 30 


+ ((A y ) 2 ) 3 2 fv 2 t ((A 0 ) 2 ) 3 2 f sin 0 
2 v 2 3 v 2 2 sin 9 30 2 


( 3 ) 


Correspondence Between the Random-Walk and Fokker-Planck Equations 

The random-walk expression (eq. (2)) describes the average results of a large num- 
ber of walks, or a large number of walks gives a numerical solution to the random -walk 
equation. An attempt will now be made to write equation (2) in a form in close agree- 
ment with equation (3) to see whether random walks can provide solutions to the Fokker- 
Planck equation. 

In equation (2) let N = 2, = v, ^ and At = Vv (where 

v = i/(v, 9) = collision frequency). By equation (1) Pj + qj = 1/2 and P 2 + ^ = */ 2, Let 


10 



2(p x - q t ) = 


2(p 2 - q 2 ) - 


\a9-1 

_ L v 


1 9( (Av) 2 ) 

+ - (A0) 

0V 

2 


1 (Av) 2 

1 9<(A0) 2 ) 

A9 Av 

V dO 

V 

\ 

1 (Ad) 2 

li 

H 

< 

^ (Av) 2 

A4 2 = -i 

1 (A9) 2 


(A cpY 


(A cp) 


5 


and 

u = fv 2 sin 0 

The random-walk equation then becomes 




/sr \ i 

U/„ ’ 2 


<Av> - 3 <(^ 2 -> + I((a 9) 2 ) + LSHl!f ((Atp) 2)) f__l 

2 2 ' \((Av)^) 


a((Av) 2 ) + 1 dy\ 


av 


2\ dv v av/ 


+ ({££) 9(( a Q) 2 > . <A9Av> [ sin 9 cos 9 /( Ay )2\| / 1 .... S((A9) 2 ) + 1 jhA 

\ 39 V 2 / \<( A g) 2 ) 39 v dej 

= .i.iU({iT) - 8< (A ^- 2) - +1 ((A 9) 2 ) + y gi- n l 9 <(A<p) 2 )\1 

y 2 3 vL \ 3v 2 2 /J 

sin ^(AS) _ 9((Ae) 2 > (A9 Av) + sin_9cos9 <( ^ )2> \1 + j <( Av ),!> jVf + j <(A9) 2 ) 8 2 f sin_9 

sin 9 89[ \ 89 v 2 /J 2 2 , 2 2 sin 9 a g2 


which differs from the Fokker- Planck equation by the terms involved with f/2 on the 
left side of the equal sign, and by the absence of the mixed derivative term on the right 
side. 

To study these equations further and to perform random walks, one must specify 
the initial and/or boundary conditions. To evaluate the Fokker-Planck coefficients, one 
must specify the type of encounter. The differences between equations (3) and (5) will be 
appraised for specific conditions. 


Physical Application 

The preceding results will be applied to magnetic mirror systems used to confine 
ensembles of electrically charged particles. The Fokker-Planck coefficients will thus 
be based on Coulomb encounters. The random walks will determine the ion distribu- 
tions inside and the end loss from a magnetic mirror system. Such a system is sketched 
in figure 3 and discussed, along with simplifying assumptions, in appendix C. The ve- 
locity space of prime interest can be best described in spherical coordinates (fig. 4). 

The polar axis is alined with the magnetic field about which there is axial symmetry. 

This B field is assumed to be uniform over the central portion of the physical space. 
The increased magnitude of B at the mirrors enters the problem by describing a loss- 
cone boundary condition in velocity space. The random walks terminate at this boundary 
as if it were an absorbing wall. These assumptions result in uniformity of f in cp and 
in coordinate space and, thus, reduce a problem in six- dimensional phase space to two 
dimensions in velocity space. 


Mirror magnet 




50 00000000000500000??03o6i0000)'300 0 0000000000 
■OOOOOOOOOOOOOQCOOOOOOOgOOOOOOOOOJCOOOOOOOOOQ^ 


Uniform field 



CD-10573-25 


Figure 3. - Magnetic mirror system. 
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The usual simplifying assumptions applied to the analysis of magnetic mirror sys- 
tems reduce the Boltzmann equation to 



( 6 ) 


for steady state (appendix C). The symbol s represents a sink (or source) density in 
velocity space. For a steady state the rate of particles being injected into the system 
equals the rate of particles entering the loss cone. Loss rates from the system can 
therefore be expressed as 


ri = - J s(v)dv 


(7) 


where integration is over the velocity space of the injection distribution. The injection 
distribution provides the initial condition for the random walk. 

To correct the random-walk results for the difference in the expressions on the 
left sides of equations (3) and (5), the following equation can be used: 


[ 
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-A 


_ =n 'f l e 2 ((Av) 2 ) 3 2 <(Ae) 2 ) 

corrected “random walk / 9 | 9 o 

z 1 av^ ae z 


<av> -. ?<(- Av > K um 2 ) 

l dv 2 


+ v_Sin 2 9< (A<p) 2 V 1 j<(Av) 2 > + i^ 

2 / \ ( (Av) 2 > 3v W 

(A e) 9((A 9 ) 2 ) <Afl Av > , sin 9 cos e ( (A<p) 2 )j ( 1 9<( a6I ) 2 > + 1 jM 

< 39 v 2 /\((A0) 2 ) 35 17 3 7. 


( 8 ) 


>dv 


The magnitude of this integral was negligible for cases of most interest. (This will 
be discussed later. ) When the mixed-derivative term of the Fokker-Planck equation is 
negligible, such as for near isotropic field distributions, n corrected is equal to that 
determined from the Fokker-Planck equation. 


Fokker-Planck Coefficients 

Expressions for the Fokker-Planck coefficients for the preceding physical applica- 
tion will now be supplied so that the step sizes and probabilities of equations (4) can be 
determined. 

The general integral expressions for these coefficients are derived in appendix D 
and can be expressed as 

/ V - V. cos 0 
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4rp 
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where (3 is the ratio of maximum to minimum impact parameter and 


r = ( Ze ) 4 ln i 3 = 0- 24 Z 4 In 

.22 . 2 

47re Q m A 


is the interaction parameter. 

In the usual problems of interest the field distribution f(v^) is not known. Also, a 
test particle is usually just a background particle with a special label. Tracing the 
paths of a large enough sample should determine test-particle distributions, which for 
a self-consistent solution should be the same as the field-particle distribution. Such a 
solution can be approached by an iteration procedure. A background distribution is first 
assumed; the coefficients are calculated; and the paths of a representative number of 
test particles are traced to determine their distribution. The distribution of test par- 
ticles in the l ^ step of this iteration are used for the field-particle distribution in the 
l + 1 step. This process is continued until the field- and test-particle distributions 
agree within the desired accuracy. The accuracy, or amount of scatter of the final data, 
in the example worked herein, was mainly set by the number of test particles used. 

In the first step of the iteration process, use of a Maxwellian distribution is con- 
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sistent with the theory of references 5 and 17. It is argued in these references that loss 
rates should be insensitive to the assumed field-particle distributions. For purposes 
herein, this may be considered a first-order solution. The insensitivity to field distri- 
bution may indicate convergence of the iteration process in a small number of steps. 

Maxwellian dis tribution of field part i cles . - A Maxwellian field distribution should 
provide an especially good first approximation to plasmas which are in near equilibrium 
conditions. An example of such could be inside mirror systems having very high mirror 
ratios so that field defects due to loss cones would be small. 

Since a Maxwellian distribution is isotropic, all first moments of 9 are zero. (An 
isotropic distribution corresponds to the diffusion approximation of ref. 5 and to the use 
of effective Rosenbluth potentials in ref. 11.) Thus (A 9 Av) is zero, and the mixed de- 
rivative term is absent from equation (3). Use of 

/ m \3/2 9 ■ mv b /2kT b 

flvpdvjj = v£e sin S b d(> b <V„ dv„ (10) 

and integrating between the limits 

0 < v^ < 00 

0 < 0^ < 7T 

0 < (p^ ^ 2n 

reduces the expressions for the Fokker-Planck coefficients to the following (see appen- 
dix D): 
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<(A<p) 2 > =-L- <(A0) 2 ) 
sin 2 0 


(lid) 


4l ^ ^ (j»X vr - l e-” V /2kTb . (l ^)erf 

9 In (5 \2KrJ »m v V mv 2 / 


Using equations (11) in (3) finally reduces the Fokker-Planck equation to 


(He) 
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An interesting relation among the coefficients of equation (11) is 


— — (v 2 {(Av) 2 )) - v<(A9) 2 ) = - (Av) 

9 / 9 


Use of equations (6) and (13) reduces (12) to equation (18. 7) of reference 5. 

Equation (12) is especially useful because there is an approximate analytical solu- 
tion to it in reference 5. The random-walk results will later be compared numerically 
to this solution. 

Using equations (11) reduces the random-walk equation (5) to 


m .iLv) . a<(Az)!) +V ((A0) 2 ))/— L — 

W c 2 |V av /\<(Av) 2 > V 3v/ 

- iiU( <AV > -iiiM^ + V ((A0) 2 ))l -<(M!i±(fcos0) + iiM!>i!fx!. + ii^ (i4) 
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In terms of dimensionless variables the average value per collision of the variables 
in equation (11) can be conveniently written as 


AV = f hi ff/ T b 


^o 2 0 2 
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and 


v Q is some reference test-particle velocity. The quantity is always refer- 


enced to and velocity v is referenced to v Q ; thus, v is actually referenced to 
through E^. If E^/T^ = 3/2, then by the kinetic definition of temperature E Q 
(p. 37 of ref. 18) is the average field-particle energy, and V is the ratio of the test- 
particle velocity to the rms field-particle velocity. Thus for most applications 
herein E^/T^ was set equal to 3/2. When studying the case of injection of a mono- 
energetic beam of test particles into a background ensemble, however, it was conven- 
ient to set E q equal to the beam-particle energy. Thus for example, when the beam 
particles had an energy of 10 times the average field-particle energy, E 0 /T^ was set 
equal to 15, making the initial value o f V eq ual to 1. 

Plots of generalized step sizes V (AV) 2 /3^ and vV(ag) 2 

yin /3 j for a Maxwellian field distribution are shown in figure 5. These curves were ob- 
tained by use of equations (15a) to (15c). 

Co mpu ting time and impac t pa rameter cu toff distance . - Setting the maximum im- 
pact parameter at the Debye length and the minimum value at a distance to cause a 90° 
deflection results in 


J3 


Debye 


0. 49x10 18 Te, 

Z 2 ^ 


3/2 


(16) 
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-3 

for (in Kev) and n^ (in m” ). Numerical values are given in table 5. 1 of refer- 
ence 12. 

If the maximum impact parameter is cutoff at the Debye leng th as in e quation (16), 
the average deflection which a particle undergoes per collision ^(A9)^/ v is extremely 
small, as can be deduced from figure 5. Thus the corresponding computing time re- 
quired to trace representative samples of test particles would be prohibitively large. It 
was observed from the elementary one-dimensional (and nonbiased) random walks of 
reference 8 and also from preliminary computer results that the average number of 
steps to reach a distance m steps away from a starting point was nearly proportional 
to m . The time to take such a walk would thus be nearly proportional to m /v- For 
a random walk in the 9 direction the distance reached per unit time would correspond 
to y(A0) , which is simply the diffusion coefficient ((A0) ). 

Solutions to processes described by the Fokker- Planck equation are dependent on 
impact parameter cutoff distance only through the Fokker- Planck coefficients. As shown 
by equations (9) and (11), these coefficients are only weakly dependent on impact param- 
eter ratio /3. It thus appears that a particle taking the fewer number of larger steps 
determined by a lower /3 cutoff would reach given locations in velocity space in about 
the same time as if it took the larger number of smaller steps determined by a higher /3 
cutoff point. The final walk time (loss rate) would possible be weighted by a function of 

fi- 

As part of the following computer study, the effect of £ on loss rates from and 
distributions inside magnetic mirror systems was determined. Cases were selected 
that could be solved analytically for comparison with the random-walk results. 

Numerical Procedure 

The procedure by which the random walks were simulated by the computer and the 
manner in which velocity distributions inside and end losses from the magnetic mirror 
system were obtained from the walks will now be described. 

A computer flow diagram and a representative FORTRAN program are presented 
in appendix E. Either the Lewis IBM 7044-7094 or 7040-7094 direct coupled systems 
were used to perform the calculations. 

Initial conditions for each random walk were determined by selecting at random 
from prescribed injection distributions. The walks were terminated at the loss-cone 
boundaries 9 and tt - 9 where 

L 
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( 17 ) 


0 C = sin- 1 


B/B 
c' o 


and B c /B q is the ratio of the magnetic field at the mirror to that in the central uniform 
region (ref. 12). 

In any of the steady- state models selected for study herein, the average number of 
particles injected per unit time equals the rate at which particles enter the loss cone. 
The walks were studied one at a time, and a large enough number of them were traced 
to form a representative statistical sample. 

The velocity space to which the test particles had access was divided into zones . 

F or each of these zones values of probabilities p^ and pg, step sizes V«Av) 2 > / v and 
^((A 9)2) / V) and collision frequency were computed for the field distribution of inter- 
est. These quantities were stored in matrix form and called for as the test particle 
progressed about the velocity field, as described in appendix E. 

During each step of each walk, a random number was selected for each of the two 
components (V and 9) and compared with probabilities p^ and pg, respectively, to 
determine step direction. These random numbers were selected from a set of numbers 
uniformly distributed over the interval from 0 to 1 (ref. 10). Whether the V compo- 
nent of the step were in the positive or negative direction depended on whether the first 
random number was greater or less than p^. In like manner, the second random num- 
ber determined the direction along 9. 

Because step sizes were so small, the test-particle locations were not recorded 
every step. They were, instead, tallied at constant time increments At- Each time 
increment included Ati >(V,9) steps. A count was kept of the number of times a test 
particle was found in each of the velocity zones (boxes). A large enough number of 

test particles A were walked to determine a velocity distribution. The time incre- 

e max 

ments were selected so that usually no less than 10 nor more than 100 steps were taken 
between tallys. The number of At's required to reach the loss cone was recorded to 
determine time of confinement and thus loss rate (appendix E). 

After each set of A max walks, new step size and probability matrices were deter- 
mined with the test-particle distribution used as the field distribution. The elements of 
the tally matrix were reset to zero and the process repeated, if desired, for the l + 1 
step of the iteration. 

Numerous computations were made to explore the parameters of impact distance 
cutoff /3, sample size A max , and iteration index for suitable field convergence 
l = ^ max - External conditions such as mirror ratio and injection distribution were 
selected either for reduction of computer time or for comparison of random-walk re- 
sults with cases having analytical solutions. 
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RESULTS AND DISCUSSION 


Influence of Impact Parameter 


The scattering angles, thus step sizes, descriptive of Coulomb encounters are, in 
general, so small that even to evaluate the effect of /3 with a reasonable amount of 
computing time, required selection of very short walks. Short walks were obtained by 
using low mirror ratios (B C /B Q » 1). The initial condition for this study was set to 
simulate injection of monoenergetic particles normal to the magnetic field 6 = v/2. 
Computations were carried out only for 1 = 1 and a Maxwellian field distribution. 

A generalized parameter descriptive of the rate at which particles first reach pre- 
scribed conical boundaries in velocity space is plotted in figure 6. These boundaries 
(0 and n - 9 c ) were selected so that no set of (either 100 or 1000) particle walks took 
more than 15 minutes of computing time. At least three (3 values were used for each 
boundary condition. The mirror ratios (eq. (17)) became extremely close to one as 
b max approached values representative of Debye lengths; in fact so low that these re- 
sults have little direct application to mirror systems of interest. These results may be 
viewed as determining the average length of walk required to first reach a certain conical 
boundary in a Maxwellian field. 

The regions of velocity space covered by the walks were so small that the step sizes 
and probabilities could be assumed constant over the length of the walks. This, in turn, 
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Figure 6. - Effect of impact-parameter cut-off distance on rate of reach- 
ing prescribed 0 distance away from starting point at ttI2. Singly 
charged particles of mass number 2; ratio of reference energy to 
field temperature, 15. 
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permitted an analytical solution of the Fokker-Planck equation (appendix F) which is 
shown by the line on figure 6. 

This analysis gave the result that the loss -rate parameter 



log j8 


was equal to a constant. This was verified over a range of /3 (from lO'^ to 10 
by the computer results to within a spread of loss parameter of about 10 percent. 

Some typical test-particle distributions are shown in figure 7. The marginal dis- 
tribution in V 


F v (V) ■ 


vM 


n. 



i(9, v)sin 9 d9 = 



F (9, V)d 9 


is a delta function in the analysis. It is sharply peaked in the computer results, nearly 
all of the points plotted lying inside 0. 998 < v/v Q < 1. 002. The marginal distribution in 
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F n{9) 


sin 9 




f(0, v)v 2 


dv = 



F(0, V)dV 


for both the theory (appendix F) and the computer experiment can be represented by 
straight lines and agree well with each other. 

The theoretical distribution from appendix F is 
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(a) Confined distribution. 



Figure 7. - Test -particle distribution for short walks and constant coefficients. Ratio of reference energy to field 
temperature, 15; sampling size, 1000; impact-parameter ratio, 10 4 ; distance between start and finish of walk, 
0. 00208 radian. 
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The theory and computer experiment appear to be in satisfactory agreement for pur- 
poses herein and encouraged applying the random-walk method to more complicated and 
practical models. 


Comparison With the First-Order Solution of Budker 

Analyt ical solution . - Among the very few analytical solutions of the Fokker- Planck 
equation applied to a magnetic mirror system is that of Budker (see refs. 5 and 19). 

This analysis makes the usual simplifying assumptions of appendix B plus the assumption 
of separability of the test-particle distribution functions; that is, it is assumed that 

f(0,v) = f 0 (0)f y (v) (19) 

inside the mirror system, and that the injection distribution (source) function is 

s (0,v) = s 0 (0)s v (v) 

-mv 2 /2kT, 

Solutions are then sought where f y (v) is of the form e . For injection of 

particles perpendicular to the B field the solution of reference 5 is 


i(9, v) = Ce 


-mv^/2kT^ 



-mv^/2kT. 9 

s(0, v) = Ce D <(A0r>cos0 6 


H) 


( 20 ) 


( 21 ) 


where the constant C can be used for normalizing. Results were reduced to this form 
by approximation of an integration over 9 by the mean value theorem. The end-loss 
rate becomes 


n = §1 ( Ze ) 4n b ln <6 3 in(V2 + l) - Vi 

(47T6 ) 2 (kT) 3/2 (Ll_\ 
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cos 9, 
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or 


nT 3//2 Vm _ 3 72xl0" 31 
Z 4 n 2 log (i 

In this solution the Fokker- Planck coefficients were based on a Maxwellian field 
distribution. It is thus like a first-order solution. 

No energy gain or loss from the system was considered other than that added by the 
injection process or carried away with the escape of particles through the mirros. Thus 
the mean energy of the particles injected E in j should equal the mean energy of the par- 
ticles escaping for a steady-state solution. Integration of the product of equation (21) 
times 1/2 mv , and then division by ri gives the mean energy per test particle to be 

E inj = 43 ^ T b 

Thus the test particles are injected at a low energy to fill in the deficiency in velocity 
space due to the predominant escape of low-energy ions. 

Random -walk solu tion. - Some preliminary insight into the random walk can be 
gained from inspection of the step sizes and probabilities of taking steps in the various 
directions. Step-size parameters for the V and 0 components and for a Maxwellian 
field distribution are plotted in figure 5. Values ur^d in a typical set of walks are given 
in a matrix form in the computer listings presented in appendix E. 

For a given field distribution, path lines, as shown schematically on figure 1(b), 

can be envisioned. Use of the step sizes V<(av) 2 )/v and V<(A0) 2 ) /v from figure 5 
permit study of the path lines for a Maxwellian field distribution. The magnitude of the 
slopes of such path lines increases with an increase of V because of the pronounced 

decrease of V(A0) 2 with V. It is apparent that the path line pattern is most conducive 
to the escape of particles at low V. It is in turn relatively easy for a particle in the high 
slope region to go to still higher velocities where escape is more difficult 

The bias terms are equally important in the test-particle behavior. Entering into 
the bias term (p^ - q^) in the V direction (eq. (4)) is a dynamic friction term AV. 

This term (discussed in ref. 7) is negative and tends to retard the test particle to zero 
mean velocity. Opposed to influence of AV are the positive valued expressions 
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X (a9 )2 + (Atp)^ - — a « AV ) 2 > 

2 2 v 3V 

which tend to diffuse the particles to the thinly populated regions of space at high V. In 

o 

spherical coordinates the available volume in velocity space varies as V , making a 
large amount of room for particles at high V. 

Similar terms appear in the (pg - q2^ which occurs in the 9 direction. All 

terms except the last (sin 9 cos 0/2)(A<p) 2 , are zero for a Maxwellian field distribution. 
For higher order solutions ( l > 1) the A0 and -A 9 AV/V terms tend to move test par- 
ticles toward 9 = 7t/2. These terms thus help confine the particles. Opposed to this 

are the (A 0) 2 } /39) and [(sin 9 cos 9)/ 2] (A <p) 2 terms that make it easier for the 

particle to escape. 

The initial locations of the test particles were determined by selection of random 

numbers from a distribution to simulate the injection velocity pattern of equation (20) as 
shown in appendix G. The walks were always started at the 9 = n/2. 

Preliminary computer results verified that the test particles which were injected 
at low energies (according to eq. (21)) had a good chance of escaping with a relatively 
low number of steps. If, however, they remained inside the mirror system and by 
chance diffused to the high velocity region, their rate of escape was much reduced, even 
though the collision frequency was increased at high V. The test particles rush to such 
higher velocity zones at the expense of energy given up by the field particles. The field 
temperature was thus corrected for the energy exchanged with the appropriate test par- 
ticle after each walk (appendix E). 

Com parison of results. - Test-particle marginal distributions from the random- 
walk procedure are compared with the analytical results of reference 5 in figure 8. In 
general, the test particles distributed inside the mirror system in much the same pat- 
tern as predicted by reference 5 (figs. 8(a) and (b)) for the first-order results ( l = 1). 

As the iteration process was continued, the test particles distributed in a slightly dif- 
ferent pattern. The peak of the distribution shifted to a slightly higher V. Beyond a l 
of 2, however, changes were within a 10- to 20-percent scatter of the data for the usual 
sample size A max of 500 particles. No attempt was made to optimize sample size 
against /3. The parameter /3 was usually selected for completion of A times 
l max walks in 20 minutes of computer time. 

Starting the 1 = 1 calculation with the field-particle distribution in the 9 direction 
cutoff at 9 and it - 9 appeared to give results as good as starting with a full Max- 
wellian distribution and completing the 1=2 step. By so doing the final number of 
steps in the iteration process (^ max ) could be reduced by one. 

The assumption of separability (eq. (19)) of the distribution function is often used to 
simplify analytical approaches (refs. 3, 5, and 19). This appears justified for the typi- 
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Generalized marginal distribution in 6,(1 - 28 c /7riF 0 (6) Marginal distribution in V, F v (V) 



Iteration 

Field distribution 

Mirror ratio, 


number, 

l 

in 0 direction 

B c /B o 

o 

1 

Full Maxwellian 

1.5 

o 

1 

Cutoff Maxwellian 

1.2 

□ 

2 

Cutoff Maxwellian 

1.2 

A 

3 

Cutoff Maxwellian 

1.2 




Figure 8. - Comparison of random-walk results with analysis of reference 5. Ratio of reference 
energy to field temperature, 1.5; sampling size, 500; atomic number, 1 ; impact -parameter 



(d) Distribution of particles as they first enter loss cone. 
Figure 8. - Concluded. 


cal set of data on figure 9. The open data points were obtained from the test-particle 
distributions of the random walk. These are a normalization of the KTJK(J, K) matrix 
(appendix E). The solid symbols denote the product of the random -walk marginal dis- 
tributions in 9 multiplied by the random -walk distributions in V. The lines are the 
analytical results of reference 5. The shift of the peaks of the computer results to values 
slightly higher than the analytical results is not due to the separability assumption. Good 
agreement is shown between the two studies for the outlet (loss cone) as well as for the 
inlet or injection distribution (figs. 8(c) and (d)). Thus the injection distributions of 
reference 5 appear quite satisfactory for a steady-state solution. 

A comparison of the random-walk loss rates with the results of reference 5 is 
shown in figure 10. Agreement is close at a mirror ratio of 1. 2. At higher mirror 
ratios, however, the loss-cone boundaries are extended enough that the chance of a par- 
ticle walking to high V fields before escape is considerably improved, accounting in a 
large part for the reduction of loss rate with increase of B c /B Q . 

The rather steep decrease of loss rate with mirror ratio obtained from the random 
walks quite closely follows the l/log(B c /B Q ) trend predicted in references 3 and 13. 

The results of reference 5 differs by the inclusion of cos 0 C in equation (22). 
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F(0, V) or F n (0)Fw(V) 




o 


Random walk with full Maxwellian in 
direction 



3.3 2.0 1.5 1.2 

Mirror ratio 

Figure 10. - Comparison of random-walk loss rates with the ap- 
proximate analytical solution of reference 5. 

The correction to the source term consisting of higher order terms in equa- 
tion (8) was evaluated numerically for a Maxwellian field distribution. For this case, 
equation (8) reduces to 


n corrected n random' walk 


0. 003x10" 30 z4r2 1o S ft 
m l/2 T 3/2 


This correction, independent of mirror ratio is negligible in the random walk results on 
figure 11. 

Computer time expenditure. - Each complete set of A walks and l steps in 
the iterative procedure required less than 20 minutes of computer execution time. It 
appears that this random-walk procedure could be extended to include special or time 
coordinates, or both, and work more general problems with reasonable computer time 
expenditure. 


CONCLUSIONS 

From a study of computer simulation of random walks of charged particles through 
ensembles of field particles inside magnetic mirror systems, and by comparison with 
results obtained from the Fokker-Planck equation for cases that could be solved analyt- 
ically, the following conclusions were reached: 

1. Computer time could be greatly reduced with no noticable reduction in accuracy 
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by a selective sampling of random numbers from the low impact-parameter range of 
the binary collisions used herein. 

2. Test-particle distributions were insensitive to the assumed field-particle dis- 
tributions. 

3. The assumption of separability of the velocity distribution appears suitable for 
application to mirror system analyses. 

4. End losses determined by the random-walk method were in close agreement with 
an approximate analytical solution of the Fokker-Planck equation at low mirror ratios. 
At higher mirror ratios the loss-cone boundaries increase permitting more test par- 
ticles to reach the higher regions of velocity space where escape is much more difficult. 
The loss rate varies inversely as the logarithm of the reciprocal of the mirror ratio. 

5. Each complete calculation, for a given set of initial and boundary conditions (in- 
jection velocity distribution and mirror ratio), was completed within 20 minutes of com- 
puter time. The study was limited to steady-state problems in velocity space. It ap- 
pears that this random-walk procedure could be extended to include spatial or time co- 
ordinates (or both) with reasonable computer time. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, October 20, 1969, 

120-27. 
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APPENDIX A 


SYMBOLS 


[The International System of Units (SI) will be used throughout with the exception of 
temperature T, reported in keV, and the corresponding Boltzmann constant k, re- 
ported in J/keV. ] 


A 

a 

a. 


a 

B 

b 

C 

E 

E 


mass number 

determinant of elements in metric tensor 
elements of metric tensor 

elements of cofactor of metric tensor 

magnetic field 

impact parameter 

normalization constant 

energy 

mv^/2k 

electric unit of charge 


v v sin 9 

F(0, V) normalized probability density times scale factors, f (v, 9) 


FJ9) 


normalized marginal distribution in 9, J F(9, V)dv 

J o 


r 


Fy(V) normalized marginal distribution in V 


I 


ir-9. 


F (9 , V)d9 


g 

H 


k 

l 

m 


probability density 

Rosenbluth potential = / f ^ u d? b 
Heaviside unit function 

Rosenbluth potential = - dv^ 

Boltzmann constant = 1. 6x10" J/keV 
step number in iteration process 
mass 
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N 

n 

n 

Pi 

«i 

q u 

S HV 

s 

r,R 

T 

rp/i 


number of dimensions 
number density 

number of particles lost per unit time per unit volume of coordi- 
nate space 

th 

probability of 1 step in the positive direction along the i coordi- 
nate 

th 

probability of 1 step in the negative direction along the i L coor- 
dinate 

general curvilinear coordinate 

r -1 <Aq M Aq y ) 

number of particles injected per unit time per unit volume of phase 
space 

random numbers 

temperature 

r -1 <Aq M > 


t time 

u magnitude of relative velocity between a field and test particle 

u(£p | 2 , . . | N ) probability of particle being located at grid location 

*1» $2’ * ‘ '» % 

Y v/v Q 

Y random number 

v velocity 

Y some arbitrary collisionally dependent quality 

Z atomic number 

/3 impact-parameter ratio, b ni ay /b Tn - n 


(Z e ) 4 In p 

r interaction parameter, = 

,22 

4ne 0 m 

5 delta function, also used to denote finite difference 

e azimuthal angle between orbital and fundamental collision planes 

e capacitivity of vacuum 
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Tj number of components changed per encounter 

d angular distance from polar axis 

$ r angle between test and field-particle velocity vectors 

v collision frequency 

| arbitrary coordinate in phase space 

a differential scattering cross section 

cp azimuthal coordinate 

Q, solid scattering angle 

Subscripts: 

b field particles 

c collisional value, also used to denote loss cone or value at mirror location 

e electron 

inj inj ection 

o reference value 

Special symbols: 

— average per collision or per step 

<) time average 



APPENDIX B 


TRANSFORMATION OF THE FOKKER-PLANCK EQUATION FROM 
RECTANGULAR CARTESIAN TO SPHERICAL COORDINATES 

The Fokker- Planck equation can be written in rectangular Cartesian coordinates as 

2 

(■£■) = -_L (f<Av^» +i_J (f<Av^ Av 1 ')) 

\9t/ c ay M 2 9v m dv v 

Superscripts are used in this appendix to denote component direction, leaving subscripts 
to denote covariant differentiation. Summation convention of repeated indices will be 
used. The subscript c on 9f/9t means changes of f due to collisions. 

Following the procedure of appendix IV of reference 11, this equation can be written 
in covariant form as 


r'Y — ) = -(fT M ) + — (fS^ v ) 

\9t/ c ’ M 2 ’ 

12 3 

which is valid for any set of curvilinear coordinates, q , q , q . 
The quantity r is an interaction parameter, 


T M = a ^(h j = < Aq^> r _1 


(Bl) 


and 


S l* v = a MW a yT (g ) = <Aq^ Aq v > r" 1 (B2) 

, Ct)T 

The terms TS^ V are dispersion coefficients, and TT^ is related to dynamic fric- 
tion. 

12 3 

For spherical coordinates q = v, q =9, and q = <p. The elements of the 
metric tensor are 

a n = i 
a 22 - v2 
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2 2 

a 3 g = v sin 0 


a. = 0 

liv 


if n * v 


and the elements of its conjugate are 


a 11 = 1 


22 1 
a = — 


.33 


v 2 sin 2 0 


3 .^ = 0 


if [1 * v 


By the procedure of covariant differentiation 


( fx M) 




where the Christoffel symbol j ^ 1 with repeated indices equals and a is 

Va gqff 

the determinant of the matrix of elements a . For spherical coordinates 


4 2 

a = v sin 0 


Thus, 


r(fT M ) = — JL (v 2 f( Av) ) +— — — (f sin 0(A0>) + JL (f<A cp)) 
’ ^ y 2 3v ' sin 0 30 3<p 


In like manner 


rttsn - 1 3 2 (yifS^^) , 1 3 
’ V* Sq^ dq v V* 3q y 
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and for spherical coordinates 


-(vifcK") - 

dq v 



vv/3?l 

Sine fS^M— ^ 
2 


9 q 


M 


3a 3a \ 
+ vv _ j 

3q w 3q" ' 


A [v 2 sin 9 f (vS 22 + v sin 2 eS 22 )l + — [sin 9 f(2vS* 2 - v 2 sin 6 cos 9 S 22 ' 
3v l j 30 l 


+ [2fv(sin 9 S* 2 + v cos 9 S 22 )] 


where use is made of the symmetry of S^. Finally the Fokker-Planck equation in 
spherical coordinates becomes 



- — — (v 2 f <Av>) i L (f Sin 9<Ae)) - — (f<A <p)) + J_ 9 2 v 2 f< (Av)^) + 1 3 2 f sin ^(Ag) 2 ) 


.2 3v 


sin 9 d9 


d<p 


2v 2 3v 2 


2 sin 9 ^^2 


2 2 2 

+ “ ) + __1 d ^ y 2 g j n A 9)) + — (v 2 f(Av A <p)) + — - (sin 9 f(A 9 A cp)) 

y 2 3v 3 (p sin 9 39 3<p 


2 3 <p 2 v 2 sin 0 8va0 


— [v 2 f(((A0) 2 ) + sin 2 0((A<p) 2 )Y| + - —If sin 0(2(A0 Av) - v sin 9 cos 0((A<p) 2 )Y| , . 

!v 2 3 v l v /j 2v sin 9 30 L V /J (B3) 

+ i^[f(<AvA<,> + V C0S 0 <A0 A *>\1 
v 3c? L V sin 0 /J 
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APPENDIX C 


PHYSICAL MODEL OF MAGNETIC MIRROR SYSTEM 

The physical model to which this analysis is applied is the usual simplified one, 
well described in many references; such as, for example, pages 186 to 189 of refer- 
ence 5, pages 3 to 5 of reference 13, pages 2 and 3 of reference 3, and pages 7 to 9 of 
reference 17. In essence, the magnetic mirror system provides a long cylindrical re- 
gion of uniform magnetic field (fig. 3) . The mirror regions are assumed to be rela- 
tively short so that essentially all the Coulomb scattering occurs in the central re- 
gion. Only binary collisions and classical particles will be considered. 

Once a particle enters the loss cone in velocity space, no matter where it is in 
physical or coordinate space, it is assumed lost from the system. Only scattering of 
ions off other ions is considered. Thus electrons serve only as a neutralizing and 
shielding background of particles. For study herein, no macroscopic electric field will 
be accounted for except possibly as a boundary condition in the manner of reference 20. 

The radius of the confining field must obviously be much larger than the ion cyclo- 
tron radius. This can be accomplished by suitable magnitude of the B field. The 
azimuthal symmetry of the B field makes v x B • V y f = 0 as shown in appendix A-l 
of reference 13. 

These assumptions culminate in the elimination of the gradient terms, in velocity 
as well as in coordinate space, from the Boltzmann equation. The Boltzmann equation 
thus reduces to 


Sf 

at 



+ s 


where (0f/3t) c is the change of f due to collisions, usually determined by the Fokker- 
Planck equation, and s represents a sink density in velocity space, uniform in coordi- 
nate space. 
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APPENDIX D 


FOKKER-PLANCK AND RANDOM-WALK COEFFICIENTS 
IN SPHERICAL COORDINATES 

The Fokker- Planck coefficients are time averages of various first- and second- 
order increments of velocity due to collisions. The random-walk process, in turn, re- 
quires the average per collision of these same velocity increments. The two types of 
averages, denoted by the angular brackets and bar notations, respectively, differ by 
the collision frequency v- Thus only one set of coefficients and the collision frequency 
need be derived for the physical conditions of interest. 

The time rate of change of some collisionally dependent quantity Y averaged over 
all scattering angles and all velocities of the scattering (field) particles is 

(Y> = /f b (v b ) yVua(u,0)dO dv b (Dl) 

where the subscript b refers to field particles, and a is the differential cross section 
of scattering over the solid angle n. The magnitude of the relative velocity between a 
test-particle of velocity v and a field particle of velocity v b is 



where 0 r is the angle between v and v^. The distribution function of field particles 
is denoted by f b and is normalized so that 

/f(v b )dv b = ^ 

Using Rosenbluth potential (ref. 11) defined as 

h=yf ( v b )idv b 

and 

g =/f(v b )u dv b 
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At- 

and letting Y equal the i u component of Av yields (ref. 11) 

< AV 1 ) = r i « 1, 2, 3, 

av 1 

the covariant form of which is given by equation (Bl) 

<Aq^> = Ta^Ch v ) 


For Coulomb encounters 


r _ Z 4 e 4 In j8 
47 re 2 m 2 

where B is the ratio of the maximum impact parameter b^,, to minimum, b . . 

I i max ’ min 

When Y is set equal to the product Av 1 Av 3 , it follows that (ref. 11) 

. . a 2 

(Av 1 Av 3 ) = r - i, j = 1, 2, 3 

9v x 9v 3 

and the covariant form is given by equation (B2) as 

<Aq^ AqO = Fa^a^fe ) 

When Y = 1, it is apparent from equation (Dl) that the collision frequency is ob- 
tained. Using the well known scattering relation involving solid angle fi, impact pa- 
rameter b, and azimuthal angle e 


a(J2)dJ2 = b db de 


in (Dl) gives 

v - 7r(b 2 - b 2 . )g 

v \ max min/ & 

which is an invariant and independent of the coordinate system. In spherical coordi- 
nates, the coefficients of main interest are 



<av> = r — 

0V 


<A9> 

v 2 se 


<a* > =- r *L 

V 2 S in 2 e ^ 

(as A,)=r(j5L.lls) 

y 2 \90 9v v ddj 

<(Av) 2 > = rii 

,„2 


((as) 2 ) =h(^M + v m) 

v 4 U 2 9v / 

( (A (ft) 2 ) = ( + v sin 2 * + sin 6 cos S 

v 4 6in 4 »W 2 3v 39 / 

<Av A„> = — r fj^i--- 13 ^ 

v 2 sin 2 0\ 9v 9< ? v9 W 


(AS A<p) = L_Mg_. £OS_0 ag \ 

v 4 sin 2 0' 90 9< ? sin 9 dcp / 


These expressions can be reduced one step further without specifying f^ since differ- 
entiation is with respect to the coordinates of the test particles. 

Azimuthal symm etry. - For aximuthal symmetry 9/9 cp is zero and thus (A cp) , 

( Av A cp) , (A 0 A cp) are zero. The remaining expressions are reduced as follows: Using 

cos 5 r = cos 0 cos 0^ + sin 6 sin 0^ cos(cp - <p-^) 
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and differentiating under the integral sign gives 


<Av> = -2r 


f - 

J (v 2 +vj 


V - v b cos $ r 


v“ + v b " 2vv b cos 3 r) 


3/2 


f (v b )dv b 


(9a) 


<A0> =-£i- 


2 y j*' vv b Jsin 9 cos 0 b - cos 9 sin 0 b cos(<p - (p b )J ^ — 




( 


2 2 

v + v b - 2vv b cos $ r 


3/2 


(9b) 


r/*' 

( A0 Av) = — I - 

V 


r j vv b (v - cos $ r )|sin 0 cos 0^ - cos 9 sin 0^ cos (<p - ^b>] 

\3/2 

b 3 r ) 


( o 2 

I V + V b - 2Wv, cos 


(9c) 


<(av) z > = r 


/ 


(v - v b cos 3 r )' 


(v 2 + V 2 - 2w b cos ^ (v 2 + V 2 - 2w b cos 3 r ) 


3/2 


f(v b )df b 


(9d) 


((AS) 2 ) =il 


f 


y 4 / I / 0 2 \ i/2 

■ 1 (v + - 2w b cos 3 r J 


2 v^ Jsin 0 cos 0^ - cos 0 sin 0^ cos(<p - cp b )J 

/ 9 9 \ 3/2 

(v + v~ - 2vv b cos S r j 


r&b> dv b (9e) 


<(a<p) 2 > = r 


4 .2 



v 2 - vv b sin 9 h cos (<p - <p fe ) 
sin 0 


v -1 sin^fl | („2 . „2 

\ 


v 4 + v b - 2w b cos 


1/2 


f(v b )dv b (9f) 




(9g) 


where b min was set equal to the classical distance of closest approach 
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(z e)‘ 


min 


4ne 


{! w ») 


_ (Ze)^ 

6,re o kT b 


The average value of these quantities per collision, in dimensionless form, are 

«V b )d7 b 


a v = < AV > = . 2 


2 ? V E o/ 


h 


V - v b COS 3 r 


v 2 - 2VV b cos s r ) 


3/2 


/V 5 


(D4a) 


V* + V b - 2W b cos 3 r f(V b )dV b 


_ <Afl) _ _ 9 In ff / T b > 


/ 


Wt, [sin 9 cos - cos 0 sin cos(co - <p h ) _ - 

si r p - b _ f(v b )dv b 


(v 2 + V b - 2W b cos 5 r ) 


3/2 


2 V 2 /3 2 \ 0 / 


/V y2 + V 2 - 2W b cos 3 r f(V b )dV b 


VV b (V - Vjj cos 3 r ) [sin d cos 0 b - cos $ sin 9^ cos(<p - <p b )| _ _ 


(D4b) 


A 9 AV = 


. <A e AV) _ 9 In (3 / T b\ J 

/ 2 2 \ 3 ^ 2 
(V 2 + V£ - 2W b cos 3 r j 


v 4 /3 2 V 2 \ E o/ 

//v 2 + V 2 - 2W b cos 3 r f(V b )dV b 


f 

! (V - V b cos 3 r ) 2 

f(V b )dV b 

_ <(AV) 2 ) _ 9 In ^/ T b\ 2 J 1 

f 2 + V b - 2VV b cos 3 r (V 2 + V 2 - 2VV b cos Sp ) 3/2 _ 

v 4 ^2 \eJ 

/Vv 2 + V 2 - 2W b cos 3 r f(V b )dV b 



(D4c) 


(D4d) 


(A0) 


2 _ <(A8) 2 ) 


_ 9 In /? 
4 V 4 /3 2 


/T \ 2 II 3 P+ 


- 2VV, cos s 
b r 


V 2 V^|sin 9 sin - cos 9 cos 0^ cos((p - <p b )J ' 

3/2 

b " * v v b r-" * 


~ 


(v 2 *v| - 2 W b COJ 8 r J‘ 


► f(v b )dv b 


V b - 2VV b cos 3 r «v*b 


(D4e) 
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(A (p) 


ln/3 , 

\2\ r> \E„ 


2 _ (JAVT) - 9 \ q 


2 f V - VV b' 


/ 


sin 0 sin 0^ + cos^0 sin 0^ 


sin 0 


cos (<p - <p b ) 


V 2 + V 2 - 2VV b cos 3 r 


4 V 4 /3 2 sin 2 © 


V .... 

/ Vv 2 + V 2 - 2VV b cos 3 r f(V b )dV b 


-f(V b )dV b 


(D4f) 


and 


v = 


4r/r 

9 In /3 


tergr / v 


v 2 + v£ - 2VV b cos « r f (V b ) dv b (D4g) 


where 


V = 


V, 

v b=- 

v o 


E = - mv 2 = kE* /3 » 1. 


o 2 o 


(D5) 


Maxwellian distribution . - A popular simplifying assumption is that the field par- 
ticles are distributed in a Maxwellian spherical distribution. This a case of the "diffu- 
sion approximation" on page 175 of reference 5 and corresponds to the use of "effective 
Rosenbluth potentials" as on page 15 of reference 19. With this assumption 

/ rv. \3/2 o -mv b /2kT b 

f b dv b = n b(^rj v b e sine b de b d<p b dv b (10) 

For a spherical distribution, 3 r can be replaced by 0 b . Integrating between the limits 

0 < v b < « 

0 < £> b < 7T 


and 


0 < <p b < 2 tt 
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gives 


g = g(v) = n b 


W 


2kT b -mv 2 /2kT b / kT 


+ v + 


7rm 


£)-* (Vi 


The remaining potential h can be obtained from 


= v2 g = !^erf 

v v \ V 2kT H J 


Since both g and h are independent of 8 and cp, it is apparent from (Bl) and 
(B2) that 

(A9) = <A cp) = (A9 Av) = 0 

The remaining terms of interest for this field distribution are 


(Av) = 


2”b r 


V 


~ -mv /2kTv, 
m ve b - erf 


27rkT v 




(Ha) 


kT 

<(Av) 2 ) (Av) 

mv 


(iib) 


((a e)‘) = 


2 \ 


Y nm v 


-mv Z /2kT b / kT- 


+ 1 ^ erf 


mv 


m 


2kT v 


(He) 


((Acp) 2 ) = — - — < (A0) 2 ) 
sin 2 0 


(Hd) 


and 


v = 


- 1) 


9 In p \ 2kT v 


JS-fv r 
2kT b y Y nm v 


-mv /2kT K / kT- 
- e b + [1 + 


mv 




(He) 
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2 

These expressions for <Av> and <(Av) ) agree with those of reference 5 and with 
those of Chandrasekhar (pp. 73 to 75 of ref. 12). These coefficients are called paral- 
lel components since they are changes in v in a direction parallel to v. Agreement 
with these references is also obtained for the perpendicular component of diffusion 
which is 


( (AVj^) 2 > =v 2 [<(A0) 2 > + sin 2 0((A<p) 2 >] 
Remaining quantities of interest are 



(15c) 


(A <p) 2 = (A 0) 2 (15d) 

sin 2 0 


and 


v = 


0. 518X10" 31 n b Z 4 ;3 2 V 
T 3//2 m 1 ' /2 




(15e) 


where B » 1. 
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APPENDIX E 


DESCRIPTION OF COMPUTING PROCEDURE 

A flow chart representing a typical calculation procedure is given in figure 11. The 
corresponding FORTRAN program and some typical results are included at the end of this 
appendix. This is preceded by a key to the computer words. A brief description of the 
program to perform the random walk studies is given in the following discussion. 

The distribution of field particles is represented by KTJK(J, K) matrices. The 
indices J and K represent partitions in V and 9, respectively. In most of the com- 
putations this was an 11 by 13 matrix, thus dividing the velocity space of most impor- 
tance into 143 zones. The initial (LL=1) field distribution was read from data cards. 

The velocity field to which the test particles had access was also divided into zones. 
These zones were identified by subscripts M and L. Step sizes, DVSQAV(M, L) and 
DXSQAV(M, L), and probabilities, TWOP(M, L) and TWOR(M, L), were computed over 
the ranges of M and L for the field distribution specified by the KTJK matrix. Inte- 
gration of equations (9) with respect to cp was by Gaussian quadratures. Integrations 

2 

with respect to V and 9, weighted by v sin 9 f(v, 9) as in equations (9), were replaced 
by summations over the J and K indices with weighting by KTJK(J, K). 

Because step sizes are so small, the test-particle locations were only recorded 
every MMAX(M, L) steps. To obtain steady-state distributions, it is desired to locate 
the test particle after each small constant increment of time, At. The number of steps 
per constant increment of time is 


RMAX(M, L) = (RMAXl)(y(M, L )) (i 

v(V2=l, X2=ir/2) 

where RMAX is the real number counterpart of integer MMAX. The value of RMAX1 
is specified at the beginning of the program and should be judiciously selected to keep 
computing time down and yet obtain accuracy. 

An alternate procedure to the use of velocity zones and matrix descriptions when 
LL=1 was to specify a Maxwellian distribution of field particles and use equations (15). 
In this way the step sizes, probabilities, and number of steps per time increment are 
continuous functions of test-particle location and are integrated over a continuous field 
distribution. A comparison of this procedure with the matrix procedure served as a 
check on the required number of velocity zones for sufficient accuracy in the matrix 
representation. Excellent agreement was obtained by use of the 11 by 13 matrices. 

3 

Suitable values of RMAX1 varied from 10 to 10 depending primarily on 

The test particles were labeled by A (used as an integer). It was found that the 
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walks of about 500 test particles (AMAX=500) was required for a representative sample. 

The initial conditions for the random walks depended on the injection distribution 
being simulated. For the example on the flow chart, V2 was always set at 1. 0 and X2 
at tt/ 2, representative of monoenergetic injection and normal to B field. This initial 
condition was used in the short walk calculations. For simulating the injection distri- 
bution of reference 5, the procedure of appendix G was used. 

The random-walk proper comprised a small part of the programming effort. In 
the two-dimensional space of interest it required the random selection of 2 numbers, 
RNR and RNS. The call, RAND (RN1), for the first random number, RN1, must be 
preceded by SAND (RN1). This sets up addresses in the congruence type of random- 
number generator subroutine of the Lewis computer library. This is a pseudo-random 
number generator inasmuch as each time SAND (RN1) is called the same sequence of 
random numbers is started again. Comparison of RNR and RNS with TWOP and TWOR 
determines in which directions the steps are taken. Step sizes and probabilities are 
held constant over MMAX(M, L) steps. After MMAX steps are taken the new location of 
the test particle is determined and tallied. Also after MMAX steps the energy transfer 
between the field and the test particle is determined and summed to that of the preced- 
ing N groups of MMAX steps. 

If the test particle has not yet reached the loss cone, it is tallied according to the 
procedure on the lower right of the flow chart. It is simply located in one of the veloc- 
ity space zones and credited to the corresponding element of a KTJK matrix. A cer- 
tain probability exists for a particle to wander about indefinitely inside the mirror sys- 
tem without being lost. A walk length of 5000 MMAX steps was therefore set as a limit, 
after which the walk would be terminated and the next particle selected. The number 
of such walks was labeled IBAD. 

When a test particle reaches a loss-cone boundary, the flow in the chart moves to 
the left. If it is the AMAX^ particle, the tallying procedure in the lower left is fol- 
lowed. The marginal distributions in V and 9 are then determined along with a 
counting of the total number of points (KTOT) tallied in the KTJK matrix. The sum of 
the values of N when the particles first reached a loss-cone boundary is printed out 
as SUM of NA. 

Denote the collision frequency per test particle at X2 = n/2 and V2 = 1 by 
v{v/2, 1); and define QUCOR by 


QUCOR = 



V 


m 

Mass of a deuteron 


Then by use of equation (9g) 
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QUCOR = 0. 896x10" 18 /3 2 



KTJK(J, K) 
KTOT 



_ : 

- 2V • SK • cos x j[ 


for z = 1. The average number of A T time increments per walk is SUMNA divided by 
AMAX. The average number of walks per unit time multiplied by 

^T V2y^\/ m / Ma s s of a deuteron is then equal to (QUCOR • AMAX)/(SUMNA • RMAX1), 
and the loss rate parameter can be expressed as 


nT 3//2 Vin _ 0. 578x10" 13 QUCOR • AMAX 
n 2 SUMNA • RMAX1 


After each walk, the field temperature is corrected to account for the energy trans- 
ferred with the test particle. Assuming AMAX particles in the field ensemble results in 



A 

SUMNA • AMAX 


S UMN A 

y ^ (v2 2 


N=1 


- VINJ 2 ) 


(E2) 


The correction is applied through the dependence of step sizes and probabilities on 
v\/m/2kT^ as in equation (10). 

After A reaches AMAX the resulting test-particle distributions are used to calcu- 
late new step size and probability matrices, and the iteration process in LL is continued 
until LL reaches LMAX. 


FORTRAN SYMBOLS 


A 

test-particle number 

AM 

AMAX - IBAD 

AMAX 

size of test-particle sample 

BMAX 

impact parameter ratio 

BRATIO 

mirror ratio 

CK 

cos(X) 

CRELi 

cos 9 cos 0^ + sin 9 sin 0^ cos^x^ 
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cx 

CXK 

DELTAV(J, K) 

DELTAX(J, K) 

DELTHA 

DELVHH 

DELVHI 

DELVSQ 

DELXHI 

DENOM1 

DVNSQ 

DVSQAV 

DXSQAV 

EOPT 

FDV 

FDVi 

FDVSi 

FDXi 

FDXSi 

FFVi 

FMIXi 

FNDi 

FNUi 

FPSIi 

FRICH(M, L) 
FTERM(M, L) 


cos(X2) 

CX*CK 
(AV)(J, K) 

A0(J, K) 

(A0) 2 

DELVHI* TR 

increment of V between successive J indices 
(AV) 2 

increment of 9 between successive K indices 

52 KTJK(J, K)GND(J, K) 

J, K 

SUMNA 

(V2 2 - VINJ 2 ) 

N=1 



E !/ T b 

1. 21/(SIB*EXP(Z*Z)) 

Ti/Ui when i = 1, 2, 3, 4 

TT2 2 /Ui 

V* SRE Li/Ui 

1.0/FNUi - Y*SRELi**2/U2 
FNUi* FDVi* FDVSi 
V*SRELi*FDVi 

SQRT(1. O+Y-2. 0*V*SK*COS(7rx i )) 
SQRT(W+Y-2. 0*V2*V*(CRELi)) 
(Ti+V*CX*SRELi/Sx) /FNUi 
2. 25*S*NUM9/(W*NUM1) 

-6. 75*S*NUM6/NUM1 
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FTHEi 
F (THETA) 

V*SRELi*(V2-2. 0*CRELi-3. 0*V2*Y(SRELi/FNUi)**2) /Ui 
KMARGX(K) print out 

FV 

ERF(Z) -R 

F(V) 

KMARGV(J) print out 

GDV 

y]w.FDVi( X -) 

i 

GDVS 

X)w.FDVSi( Xi ) 

i 

GDX 

^ Wi FDXi( Xi ) 

i 

GDXS 

J^w i FDXSi(x i ) 

i 

GFV 

^ Wi FFVi( Xi ) 

i 

GMIX 

2^w i FMIXi( Xi ) 

i 

GND 

J]w i FNDi( Xi ) 

i 

GNU 

^w.FNUUxj) 

i 

GPSI 

JV.FPSUfxj) 

i 

GTHE 

^ FTHE i (Xj) 
i 

IBAD 

number of walks discarded due to N reaching value 5000 

IX 

number of steps in positive 8 direction 

IY 

number of steps in positive V direction 
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KJ(J) 


KMARGV(J) 

KMARGX(K) 
KTJK(J, K) 

KTOT 

LJ(J) 

LL 

MMAX 

N 

NA 

NUM1, NUM2, 

QUCOR 

R, RN1, RNR, 

RAND 

RMAX 

RMAX1 

S 

SAND 

SIB 


number of points in J th interval of injection distribution 
in V 

^KTJK(J, K) used for marginal distribution in V 
K 

F KTJK(J, K) used for marginal distribution in 9 

tally of points in J increment of V and K Ln increment 
of 9 for field distribution determination 

£ KTJK(J, K) 

j,k 

number of points in J Ln increment of V in loss-cone 
distribution 

step number of iteration process 
number of steps between tally points 
number of groups of MMAX steps 

number of groups of MMAX steps when particle reaches 
loss- cone boundary 

. . . , NUM9 used in weighting functions by KTJK and summing over 
J and K 

collision frequency at V2 = 1. 0, X2 = n/2 
RNS random numbers selected from uniform distribution 

call code for random numbers 
real MMAX 

MMAX when V2 = 1. 0 and X2 = tt/ 2 
In /3 

(E^/T b ) 2 

initial call code to set up addresses in random number 
generator 
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SK 

sin X 

SRELi 

sin 0 cos 0 b - cos 9 sin 0 b cos^x^) 


AMAX 

SUMNA 

T, NA 


A=1 

SX 

sin(X2) 

SXK 

sin(X2) sin(X) 

Ti 

V2-V*CRE Li 

TR 

1 - AT b /T b , see eq. (E2) 

TTi 

V*SQRT(1. 0-CRELi**2) 

TWOP 

2pj 

TWOR 

2p 2 

Ui 

FNUi**3 

V 

magnitude of field-particle velocity 

V2 

magnitude of test-particle velocity 

VHH 

Vffl*TR 

Vffl 

initial location on V scale for tallying 

VHT 

initial location on V scale for tallying 

VINJ 

injection velocity 

VO 

V 1. 5 T b /E^ arcsin(R 2 ) 

VOUT 

V2 value of point being tallied 

VSQFV(J) 

V 2 KMARG(J) 

VSQFVT 

Ev 2 KMARGV(J) 
J 

VSQKJ(J) 

v 2 kj(j) 

VSQLJ(J) 

v 2 lj(j) 

VSQKJT 

Ev 2 KJ(J) 

J 

VSQLJT 

Ev 2 LJ(J) 


J 


particles 

particles 
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XMIN and XMAX 
XOUT 


Subscripts: 


(V2) 2 

weight factors in Gaussian quadrature integration procedure 
abscissa location in Gaussian quadrature integration procedure 

0 b 

9 

initial location on 9 scale for tallying particles 
initial location on 9 scale for tallying particles 
loss cone boundaries on 9 
X2 values of a point being tallied 
V 2 

1.073 VO 

index on V or V2 
index on X or X2 
index on V2 
index on X2 

index on Gaussian quadrature terms 
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FORTRAN PROGRAM 


C 1 1* -NS ION KM Aj<GV'( 11) , KPARGXC 13) ,KTJK( 11, 13) ,KJ ( LI) ,LJ( 11 ) 

Cl M . NS I ON r V St Ay ( u * i 3 ]_J i;XS.Q A V (I Lt_L3 J * TWi)P.( U , 13 L. I WQR U 1 * 1 3 ) , 
L MXU< 11*13) v VSOFV(ll) 

C I MiiNS I DN UEL[ AV ( L L*_L3.) * CELT AX (11* 13 J , OPS I I LI, Id) »E TERM ( 11 , 131 
1FRICH (11,13) * CM I X ( 11, 13) 

INTEGER A,*MAX 
CALL SAND(RNl) 

C_JEAD I\ OF INITIAL FIELD C I S TR I OUT I dM _ 

REA t : (5, 50 2) ( (KT JK( J,K) ,J=1 ,11 ) ,K=1, 1 3) 

*EAf (_5t5Q3J (KNAKCVCJ ) ,J= U 11) 

PEA 0(5* 50 4) (KMARGX(K) ,K=1, 13) 

REAL 15, 50 5 IKTU.T 
502 FCRMTl 1114) 

5jii. F.Q R m AT (. 1 1 1 *3 1 .. _ .. _ 

:>04 FCR^AT( l 315) 

505 . FORMAT ( 15 ) 

UR I r F * 6 ,625 ) 

. . UP ITE(6,.62P) (K^IKTJKIJ**)* J=l # 11 ) ,KMARGX (K ) . K=1 , 1 3 ) 

ViRI I E( 6,629) (KMARGVIJ ),J = 1,11) 

_.kR I f 5-LA, 6 3 CLlfcT QT_ .. ... . 

10 PF A l . (5 , 50 0 ) ECPT , UMAX , RMAX1 , A MAX, LMAX 
..If (. LQP.T * EG.* 0 » ) GO. ID 170 

REAu (5 f 501 ) XMIN, XMAX , VFI ,XHI , DEL VH I , DEL XH I 
_ . UR I f E ( .6 j 60 C ) EUPTjJbMAX , RMAX1 , AM AX »XMIN*XM AX . 

PRAT 10=1.0/ (SI N ( X M I N ) )**2 

J? R X.TJ: 1 6 * 6 3 1 ) VJbLI.i. X H. Lf D.tLL VJz I_j D £ L X LLl_* BRA TJLU . . _ .. . _ „ 

CC30LL“1,LMAX 
T R =_1j > 0 

$ = A L OG ( LMAX ) / (EU PT*BMAX) * * 2 

C_ I TERATIj'N J)0 LOUP FOR SELF CONS I STFNT_ F I ELD .01 STRIBUTI ON 
CC2M = 1 , 13 

C CALC L L .41.10 N. ..IE. iTJtP -S.LZ £._ £Ml_PP..at: ABJLU LTY M4I&I.CE S. 

L C? 5 M= 1 , 1 1 

V2=vHI + ( FL^AT ( M )-1.5 ) *DELVHI 
X2=XHI+ ( FLOAT (L )- 1.5 ) *CELXHI 
U= V 

CX=CCS ( X 2 ) 

_S.X =dl MLX2 1 . . 

REAL N U w 1 , NUM2.NUM3, NUN 4 , NUM 5 , NUM6 , NUM7 , NUM3 • WUM9 

. MV1=0.0 ... .. 

MM2=0. 0 

MJM3 =0 . 0 ... 

MJM4=0. 0 

, Ayjds=o,o_ . _ .. 

MM 5=0.0 

mj£/=o.o.. . .. . ... ......... 

MJM=0.0 

_ A L M i=0. A 0 . . .... .... ...... 

cf.\cmi=o.o 

£ I IN T EG.R A I LGL6L Q.VJE P G L Aid. AN G L fL_AN U. VELOCITY. B.Y QUADRATURES 

CG31K=1,13 


I 


r c 3 ? J= l * 1 1 _ _ 

IF(LL-lil3, 13,14 

13 V= (0*200+ f FLOAT I J 1-1 * 5 )*C* 200 I *SQRT ( 1.5/EOPT ) 

X = 0. 241 66 1+( FLOAT ( K ) — 1.5 )*0. 241661 

GO TO ] 5 ... . ... 

14 V=VHI+( FLOAT ( J)-1.5)*CELVHI 

X = Xhf+(FI rflTCKJ-1 - 5)»rFI XH T 

15 Y = V **2 

CK=cns c X) 

SK=SIN( X) 

CXK=CX*CK 
S X K = SX* SK 

C T MTFC R ATI FI N Q VFtt A7 IHIJTHA L ANGLE BY GAUSSIAN Q UA OR A T U R FS 

CREL1=CXK+SXK*C0S(0. 1 € 34 3* 3 . 14 16 ) 

S R E L 1= SX * C K— CX *S K * CO S { 0 * 18343* 3 * 14 16 ) 

FNUI = SQRT ( W+Y-2.0*V2*V*CREL1 ) 

FNO l = SQRT C le-O+Y— 2sO*V*SK*CCS (0*18343*3*1416) ) 

L1=FNU1**3 
T1=V?-V*0RF!_1 
TT1=V*S0RT ( 1*0— CR ELI* *2 ) 

FOXl=y*SRELl/Ul 
FD V 1=T 1 /U 1 

FD VS 1=1 T 1 **27 U 1 _ 

FDX$1=1.0/FNU1-Y*SREL1**2/U1 

FFV L=FN!t j l* F n yi *FDySl „ „ 

F PS 1 1= ( T1+V*CX*SREL1/SX) /FNU1 

FM? xi=V*SREL 1*F0VI * _ _ 

FTHE1=V*SREL1*(V2-2.0*CRELI-3.0*V2*Y*(SREL1/FNUI )**2)/Ul 
CREL2»CXK*SXK*CCSC4»52553*3e 1416) - 

SREL2=SX*CK-CX*SK*COSCO. 52553*3.1416) 

FNU2=SQ RT g w +Y-2 c Q*V2*y*€ REL2 ) . . _ 

FND2=SQRT ( 1 . 0+ Y-2 . 0* V* SK *C0 S ( 0 . 52553* 3 • 14 1 6 ) ) 

L 2 = F MU 2 * * 3 .. .... . . . .. 

T 2= V2-V*C REL 2 

1 T2=V*SQRT Cl* 0—CREL2 **2 1 . _ 

FDX2=V*SREL2/U2 

F Dy 2=T2/U2 - 

FCVS2=TT2**2/U2 

FDX S2= 1 e 0/FNU2— Y*SREL2»*2/U2 

FFV2=FNU2*FDV2*FDVS2 

FP S 1 2= ( T2 + V*CX*SREL2 / SX ) /FNU2 

FMI X2=V*SREL2*F0V2 

F THF2=V»SR EL 2* (V2— 2*0»CREL2~3*0*y2»Y*(SREL2/ -FNAJ2 )** 2 )/U 2 

CREL3=CXK+SXK*C0S( 0.79667*3. 1416 ) 

5 R E L 3= S X * C K— CX * S K* CO S (0*7966 7*3*1416) 

FIMU3 = S0RT ( VK Y— 2 • 0*V2*V*CREL3 ) 

FM03=SQRT ( 1-0+ Y-2*0*V*$K«CQS (0^79667*3* 1416) ) 

L3=FNU3**3 

T 3 = V2~V*C REL3 .... _ „ 

1 T 3 = V*SCR T ( 1 . 0-C REL3 * * 2 ) 

FDX3-V*SREL3/U3 
F D V 3 =T 3/U 3 

p n y s 3= T T 3 * * 2/ U 3 ... _ _ . „ __ 

FDXS3=1.0/FiNU3-Y*SREL3**2/U3 

F F V 3 :=F Ni l 3 ♦ F U V 3 » F P V S 3 __ 

FPS I3=(T3 + V*CX*SREL3/SX ) /FNU3 
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Ill III I 


FP I X3 =V*SREL3«F0V3 

FTHt3=V*SRFL3*( V2-2. 0 *C R EL 3- 3 . C* V2*Y* ( SREL 3/ FNU3 > * *2 )/U3 

£_R,E14=CXK + SXK*CCS( 0*9.6.025*3^.14161 ... _. .... .... 

SREL4=SX*CK-CX*SK*CUS ( 0.96 02*5*3.1416) 

F NU t = S QRT (. W+_Y_-2 • 0*V2 *V*C R.EL4 1 

FN04=SQRT ( 1 .0+Y-2.0*V*SK*C0S (0.96029*3. 1416) ) 

L4 = FMU4» «3 _ .... ...... _ _ 

T4=V2-V*CREL4 

TT4 =V»SCRT ( 1 . 0-C REL4** 2 ) . . .. ... ..._ 

FDX4=V*SREL4/U4 

F0V 4=T4/U4 . ... . . 

FCVS4=TT4**2/U4 

FDXS4=1. Q/FMJ4-Y»SREL4*«2/U4 

FFV4=FMU4*FuV4*FCVS4 

FPS 14= ( T4+V*CX»SREL4/SX ) /FNU4 ...... . ._... .. _ 

FNIX4=V*SREL4*FCV4 

F TH'.4 = V* SRFL4* ( V2. - 2« O*CREL4-3«_0* V2* Y* ( SREL4/F.NU.4 1**2) /.U4 .. ,. _ 

G NU =0. 362 6 8*FNU 1 + 0.31 371 *FNU2+0. 22238*FNU3+0 . 10123*FNU4 

C NO = 0. 36268»FiM01 <-0.3 1371 *FND2+0. 22?38*FND3+0. 1 0123» FND4 .. 

GDV=0. 3626R*F0'/1+0.31371*F0V2 + 0.22238*F0V3 + 0.10123*FDV4 

G.Q.X_=9. 36268*FDXl+0. 31271 *FDX2+0_. 22238*FDX3+Q..10122*FDX4. 

GDVS=0.3626o*FDVSl+0. 3 13 71 *F CVS2 + 0 . 222 38*FDV S3 + 0 . 10123*FDVS4 

£0X3=0 . 3 6.26 <3* FOX S 1+0 • 3 13 71* F CX S2 + Q .2 223 8* FD.XSJ+.0 . 1Q123*F.DXS4 

G FV=0.3626R*FFV 1+0.3 1371 *F FV 2+0 . 22238*FF V 3+0 . 1 0 l 23*FF V4 
GPS 1 =0. 3AZhb *F PS 1 1+0 . 3 1371 *E P S 1 2 + 0 . 22238+F PS 1 3 + 0 . 1 Q.123*FP SI 4 . 
GPr<=0.3 62 68*FMIXl+0.313 71*FPIX,' + 0.222 38*FPlX3 + 0.10l2 3*FMIX4 

GTH.. =0.36268»FTHE1 +0.3 1371* FTHE2+.Q «22228*ETHE.3.+ 0* IQ123*F T.HE4 

NU W 1=NUM1+FL0AT(KTJK ( J ,K ) )*GNU 

&_U. v 2=NiJ!'2 + FL04TIKTJK 1 JiK U *.GDV _. . _ 

Ml* 3=NUM3+FLCAT( KTJK (J,K ) ) *GCX 

Ml5i / i = NUtl4 +FLCA F (KT JKLJiK.).)_*GCVS. . _.. ..... . ..... 

MM = NUM5+FIGAF ( KT JK ( J ,K ) )*GDXS 

MJ£j2=NUy6_ + FLQAT..(KT.JK( J,K)1*GFV ... ... . .. .. 

MJM / =NUM7 + FLCAF (KTJK { J,K )) *GPSI 

MJP.3 =NUP8 + F.LC A.T { K J JK LJ * K ) ) *GP I X 

MP3 =MIJP9 + FL CAT { KTJK ( J,K ) ) *GTHE 

C^L£Pl = DE.N|.CPl+fcL.CAXiKT liajjjLU.*GND 

32 CONTINUE 

31 CONTINUE .. ... . ....... 

PPAX=RMAX 1*\UM1/0EN0M 1 

CELT AV{M.L)=-4.50*S»NLP2/NUM1 

CEL fAX(P,L )=-4.50*S*NLP3/( V2*NUM1) 

C£LVSQ = 2.25-*S*:NUfrl4/NUPl .. .. . ... 

CELTHA = 2. 2 5*S*NUM5/( W*NUP1 ) 

E.Y3CAV (P,L ). = SURT(jiEL12SlLl__.. 

CXSGAVJM, L )=SORT( DEL TEA) 

F.T.ER.iyil M » L ) = -6.» 75*3 *NUP6./.NUM1 ... .. 

CPSI{M,L)=?.25*S*NUM7/(W*V2*SX*SX*NUP1) 

_£IiI3_(M ,.L i =2.»2iL* S*NUM 8 Z1V 2*!IUP 1 ) ... _ 

FR ICH( M, L )=2.25*S*NUM9/( w*NUPl ) 

_ IhCPlFULJ =0.5* (1.0 + 1 0ELTAY(.M,L)-FTERHIM,L)+0.3*V2* L0ELtHA+SX**2* 

1CPSI(M«L) ) )/DVSQAV(M,L) ) 

JUdCLL l. M ».L1.= 0 . 5.*.(J.-Q.+ 1 DELTAX.tM • L ) -FRIC.HIM » LLrCLMI XI M*L1 /Y2+.Q.5*SX*CX* 

1CPSI (M,L) ) /DXSQAV ( P» L ) ) 

(LP.AXIM.»JJ=.IFIXIRMAX+.Q^5) .. ... ... . . 
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2ft CflNTINIJF 

29 CONTINUE 

CUCnR=0-fl96F-1fi*fl MAX»«?»<;QRT!FnPT)* DENOM1 /FLOAT f K TflT \ ■■_ 

WR I TE ( 6 , 624 ) LL 

WRITF16.632) 

WRITE(6»633) ( L , ( DVSQAV (M ,L ) , M= 1 , 11 ) ,L=1 f 13 ) 

WRTTFC6.634) . 

WRITE(6,633) (L, (DXSQAV (M»L ) » M= 1 , 1 1 ) , L= 1 , 13 ) 

kR ITEC6-635) 

WRITE (6, 633) CL f !TWOP!M,L > , M= 1, 1 1 ) , L= 1 , 13 ) 

kRITFI6,636) 

WRITE!6,633) !L, (TWOR! M,L ),M= 1,11), L= 1,13) 

WRTTFI6.637) __ 

WRITE! 6,638) (L, !MMAX!M,L ) , M = 1, 1 I ) , L*l, 13 ) 

rni7J=i,ii 

KMARGV! J)=0 

K J I J ) =0 

L J ( J ) = 0 

mi 7 k=i .13 

KMARGX! K) =0 

17 KTJKCJ-KJ =0 

KTO T=0 

CVNSQ=0^0 

$UMNA=0 • 0 

1840=0 __ 

SIB=1*0/SQRT!E0PT) 

C START QF WALKING AM AX PARTICLES, ONE AT A TIME _ 

C090A=1,AMAX 

C ALL RAND C R ) ... 

C NEWTCN RAPHSON DETERMINATION OF V!R) FROM R!V), BUDKER f S DISTRIBUTION 
vc= l^225*SI£*ARSIN4R*R } J 

1 2=1. 07 3 *V O/SIB 

FV = F»F f 7 1-R 

FDV=1.21/!SIB*EXP!Z*Z) ) 

F = FV/FOV 

V2= VO-F 

IF 4ABS CF) 

2 V0=V2 

GO TO 1 

3 V I N J=V2 

VHT = 0^1 ... 

C TALLY OF INJECTION MARGINAL DISTRIBUTION IN V 

CC4J=1 - 11 

IF! VINJ.LT.VHT)GO TO 5 

_^4 VHT=yHT+Q * 15 

J = ll 

S kj ( J>=KJ ( J ) + 1 

X2=l *57079633 

K=l 

*1 VHH = VH I*TR 

CELVHH=DELVHI*TR 

C LOCATION OF INDICES IN STEP SIZE AND PROBABILITY MATRICES 

MaiFix* CV2-VHH!/DELVHH+2*Q ) 

IFIM.LT.1IG0 TO 18 

|F{M-GT=11)GQ TO 20 

GO TO 19 
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18 N = 1 

GC 10 19 

20 N = 1 i . . _ _ . . „ . _ _ _ 

19 L = I F I X { { X 2 - a F I )/CELXHI + 2.0) 

C R ANDOM *ALK PROPER „ 

I x = o 

_ i v =o __ _ 

NPX-MMAX(M,L) 

C C 2 8MM= 1 * MX „ 

CALL RANC(RNR) 

CA LL RANC(RNS) 

I F ( -’ NR • L E • T rt 0 P ( M , L ) ) G C TC 22 

I Y = I Y- 1 

C-C TO 2 3 

22 ,IY = .IY+l„ _ . 

23 IF{RNS.LE.TfcCR(M,L))GC TL 24 

I X= I X- ... 

GC TO 25 

24 I X = I X + 1 

25 CONTINUE 

C L OCATING TEST PARTICLE IN VELOCITY SPACE 

V2=V2+FLOAT{ IY)*DVSQAV(M*L) 

C VN 3Q=0_VN SQ + V2 *V 2- VI N J * V_INJ 

X2=X2+FLCAT( IX ) *DXSQAV (M ,L ) 

IF ( X2«LT.XM .1N.0*.X2«GT.XPAX)GQ 10, 86. _ 

V C U T = V 2 

X CU T =X 2 . 

X H T - X H I 

YHT=VHI ... . ... 

C FIFLC CISTRreUTICN TALLY 

CC8 LK= 1 • 13 ... „ _ „ 

IF ( XQUT.LT.XHT )G0 TO 82 

...81 > H T ■= XJdT.+ D £ LXE 1 . .. ___ 

K = 1 3 

82 CGN1INUE ._ . _ . . 

C 08 3 J= 1 t I 1 

_ IF( VCUT ♦ L T • V H T )GC TO 84. .... 

83 VH T = VHT + D ELVh I 
J=ll 

84 KTJK(J*K)=KTJK(JfK)+l 

IF( N t *GT»5C00)GG TG 80 

85 fs = N + 1 

CC TO 2 1 „ . 

86 SUMNA=SUMNA+FLOAT(N) 

C_.AC.CCLM IMG FOR FIELD ENERGY CHANGE. 

TR 1-1. 0-0 VNSC* FLOAT ( A )/( SUMN ** FLO A T ( AM AX ) ) 
IF C JR1-0. 0)67,87. 88 

87 TR-0.001 

GC IQ 89 . _ ... . 

88 TR=SQRT{TRI) 

8 9 _ . CONTINUE- . ... . „ . 

VHT^VHI 

__C.TALLV.CF PARTICLES.. FALTERING LOSS CCNE 

CC 7 J = 1 v II 

If (_y_2.LT.yHT)GE! TC 8 ... 

7 VHT=VHT+DELVhI 
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J = 1 1 . _ . 

8 LJ(J)=LJ(J)+1 

GC TO 90 

80 I0AOIBAC+1 

90 _ CCNTINUE 

91 AM=FL0AT(AMAX-IBAC) 

V$QFVT = 0,Q_ 

C TALL V FOR MARGINAL DISTRIBUTION IN V AND THETA 
CGifciJ-Ie 11 
C C 1 60K = 1 v 13 

160 KMARGV ( J) = K M A R G V (JJ + KTJK ( J«KJ 
V=VHI+( FLOAT ( J)-l .5) *CELVHI 

C CAI CU A r ION OF. HlFI D EN ERG Y 

VSCP V ( J )=V*V*FLGAT (KMARGW J ) ) 

V£G.FVT = VSQFVT + VSGFV( J ) _ . 


161 

CONTINUE 





- 

C0163K=1^13 
CC162J=lf 11 

- - - 

■ - - - 

. - -- 
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K M A RGX f K 1 = K N' A k G X ( K 1 4*K T .Ik fl.K \ 
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CONTINUE 

CE 1 6 AK= 1 -33 . - - - 





16A 

KTCT=KTOT+KMAKGX (K ) 

-kiUTEI 6,AQC1 - 
URIIE(6,A01) ( K J ( J ) , J = 1 * 1 1 ) 


-- 




UR I TF ( 6 -675 1 „ . _ . _ 

URITE(6,628) (K* { K T JK ( J , K )» J=l f 11), KM ARGX ( K ) , K=l, 13 ) 

— WR I rEt6 T 629.Ll KMARGV ( J-4-, J =1-, 1 1 ) 

kRITE!6,64l) (VSQ<=V(J ) «J=1* 11 ), VSQFVT 

.URl rEC6^63QJKTCJ „ . ~ - _ . . . _ 

URITF(6, 530)QUCCR 

U.R I TEC 6. 4 023 . 

URITE{6,A03) (LJ(J)tJ=Ifll) 

URJJEI 6.6ACJSiiMf^m _ .. .. ... 

UR I TE!6,6CC )EGPT f UMAX, RMAX1, AM AX, XMIN.XMAX 

30 CCNIINUE. . . .... . „ . . . _ 

GC TO 10 

1 7Q 5 jjnp „ _ ... 

A 00 FORMAT! 1H L , A 5X , A 1 HT A L L Y CF PCINTS IN INJECTION DISTRIBUTION 

I /30.X « 3HJ = 1 . bX t 3H J-2? bX , 3HJ-3 * 5X* 3HJ=4, 5X*3HJ=b. 5X*3HJ=6* 5X»3 HJj=Z-*_ 
25 X, 3HJ=8, 5X» 3HJ = 9,AX , A H J = 1 0 , AX , AH J = l I ) 

A Qi F CRM AT 1 lQX*5JHKJi J ) * 10X_- 1 llfi J ... . 

A 02 FORMAT! 1HL, A5X, AlHTALIV LF PCINTS IN LOSS CONE DISTRIBUTION 

1 /3QX -3HJ = 1 , bX-3HJ=7- 5X - 3KJ = 3 * 5X - 3H J = 4. 5X-3HJ=5 -_SX, 3>-iJ = 6 - r 3H ) = 7 * 
25 X , 3H J = 8 , 5 X , 3HJ=9 f AX , A H J =1 0 , AX , AH J = 1 1 ) 

AQ3 F G R v AT 1 1 0 X -_5 HL JT J J si Q X - 1 i I 8 ) 

bOO F0RMAT{F6.1,2<ElA.7) f I6,I3) 

bQI FORMAT (6CF1Q* 8} ) . .... 

600 FORMAT! 1H1, bOX, 23HRANOOM WALK CALCUL A T I ON/ 1 HO , 10X , 7HE0 '/ T =, 

L F6, 1- 5 X - 6 HE M AX = * 1PE 10* 3- 5X ? AHRMAX1 = , IPE1Q* 3-3X- SHAM AXsl+ 1 5 , b X , 

2 5HXMIN=,0PF10„7,5X,5PXMAX=,F10.7) 

62 A FORK AT tlHL , 55X* 2 OH IT ERAT 10 M -NUMBER LL = * I U .... 

625 FCRMAT!1HL,A2X, A8HT0TAL TALLY OF NUMBER OF POINTS IN EACH INTERVAL 

1 /14X* 1HK.15X,3HJ=1* 5X , 3HJ =2 , 5X* 3H J=3 , 5X * 3H J = A ? SX* 3HJ=5 , 5X s 3HJ=6+- 

2 5X,3HJ=7,5X»3HJ=8f5X,3HJ=9 f AXtAHJ=lCfAXfAHJ=lI,3Xf8HF(THfcTA) ) 

628 FORMAT f LOX t I5 ; 10X ? 1218 ) . - _ 

62 9 FORM AT { 1 IX , AHf ( V ) , 10 X , 11 18 ) 


63 



630 FCR^AT(1HL.3X .22HTCTAL ££I NLT.S .TALL LED . = ,110) _ . _ 

530 FCRi"'AT(6X»30HCCLLISION FREQUENCY PAR AKETER=, E14. 7) 

o31 FCR RAT(6X.4HVM=.F10. 7 . 5 X , 4HXH I = , F 10 . 7 . 5X . 7HD£ LVH I = . F 10. 7 . 5X . 7HDEL 

1 >H I •= »F 10. 7,5X, 13HPIRRCR RA T 1 0= , F 1 1 . 7 ) 

632 FOR MAT ( 1HL«55X»21HSTEP SIZE CVSUAV.IMjU . . _ _ ... . . . 

1 / 1 \ » 2 H L, 7x,3HM=l,8X,3Hf=2, 8X , JHM = 3 , 8X , 3H ^ = 4 , 8X f 3HW=5 , 8 X, 3HM=6, 

2 8X,3HN=7,8A, 3 HM= H,8X,3HP=9.7X.4HM=10,7X,4HM=11) 

633 FOR N AT ( LX, 12, IX, 11E1 1 .3 ) 

634 FCRPAT ( 1HL.5BX,21HSTEP SIZE CXSQAV < R! ,L ) _ . . 

1 / L a , 2 h L,7X,3HM = l,8X,3HP=2,8X,3HM = 3,8X,3HM = 4,bX,3HP = 5,8X,3HM=6, 

2 8 X,3H P = 7,fiX,3FR = 8.8X.3HP=q.7X.^HM=10.7X,4HM = il) 

635 FGRf' / AT(lHL,35X,21HPR0F.ABILITY TwOP(M,L) 

1 /IX, 2H L,7X,3HR=1 , 8 X.* 3 H 2 , 8 X ,3 HM = 3 , 8 X , 3h f* = 4_, 8_X , 3 HM = 5 , 8X, 3HM= 6. _ 

2 8X « 3HM = 7, FX,3HM = 8,8X , >HF = 9 , 7X , 4HM= io, 7X ,4HM= i 1 ) 

6 36 F C R "AT (1 HL , 5 5.X , 2 1 HPROE.A ti IL I T Y..T WOR ( P , L )_ 

1 / 1 ' ,2H L , 7a,3HM=1 ,8 X ,3HP= 2, 8X, 3WM = 3, 8X,3HM = 4, 6X, 3HM = 5,8 X, 3HM=6, 

2 8X,3HM=7,8X.3HM = 8,8X,3HP = 9,7X,4HM=10.7X.4HM = 11) 

637 FORMAT ( 1HL , 55X , 2 1 FiFR ECUEiNC Y KKAX(M,L) 

1 / 1 x ,2H L,7X ,3HK= 1 , 8.X , 3HR=2 , 8X, 3HM = 3 , 8X , 3H.M=4., 8X , 3HM=5j 8 Xj3HM=6,. _ 

2 8X,3HH=7,8X,3HM=8,8X ,3HP=9, 7X , 4HP= 1 0 , 7X , 4HM=1 1 ) 

Aid FCR'*4T(lX,.I2,iX,llIllJ . ... 

640 FCRYATt 1HL,7FSUR NA= , 1PE 14. 7 ,5X, 3HAM=» 1PE14.7) 

641 FCRRATt 1 1 X , 6h v VF. ( Y ) , 8X , 1 1F8 . ljZXutFlO.L) ... 

END 
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APPENDIX F 


ANALYTICAL SOLUTION FOR CASE OF SHORT WALKS 

Consider the walks so short that the Fokker-Planck coefficients are constant over 
the distance traveled. Assume that the walk terminates when a particle first reaches 
a prescribed 9 distance from its initial location. Effects of V are thus of second 
order and can be neglected. For a spherical field distribution, both equations (4) 
and (5) in combination with a source term as in equation (6) reduce to 

- ((A Q)^) — (sin 9 — ) = -s sin 9 (FI) 

2 d9 \ a 9/ 

If the initial test-particle location is at jt/2, then sin 9 « 1 and s = s Q 6[0 - (tt/ 2)] 
reducing equation (FI) to 

2 

- <(A0) 2 ) 11 = -S o(e (F2) 

2 d& 2 ° \ 2 ) 

This model may simulate the end loss problem for mirror ratios very close to 1. 0. 
Injection would be normal to the B field at a constant rate s . The boundary conditions 
are 


f(0 c ) = f(*r - 9 C ) = 0 


(F3) 


Using (Fl) in (F2) and integrating both sides of the equation gives 


0f _ 3f_ 
39 39 


((A0) 2 ) 



where H is the Heaviside unit function. Integrating a second time and using equation 
(F3) yields 
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m - 

dd 


(9 - 9 C ) = 0 


9=9, 


2s / s 

-VH 

/amA X */ 


<(A0)*> 


when 0„ < 9 < - 
c 2 


when — < 0 < (n - 9 ) 
2 c 




(F4) 


J 


At 9 = n - 9, 


d9 


9=9 c <(a er) 


(F5) 


Using equation (F5) in (F4) gives 


f( 0 )= 2_ ( 0 _ 0 ) 

<(A0) 2 > 


when 0 ^ 0 s — 

C 2 




— — (ti - 9 - 0 ) when - < 0 < ti < 0 


<(A0) 2 ) 

Loss rate must equal injection rate for steady state so 


r”-°, 

„= / 

«/0„ 


s (0)sin 0 d0 = s 


/>7T-0 C 

. / a( e -l)s, 


sin 0 d0 = s 


(F6) 


The number density, distribution function relation must be 


n b = 


/ 7T-0 C 


f(0)sin 0 d0 = . 


2s, 


( (A0) / 




(0 - 0 )sin 0 d0 


2s, 


<(A0) 2 > 


- — (1 - sin 0 c ) 


(F7) 
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Using equation (F6) gives 


n b ((A0) 2 > 


n = 


2(1 - sin 6 C ) 


For small — 9 


c> 


sin 9, 


=1 .(M 


(F8) 


to second order so that the loss rate becomes 


. n b <(A0) 2 ) 


n = 




(F9) 


Evaluating ((A Or) for a spherical Maxwellian field distribution (eq. (11c)) and sub- 


stituting into equation (F9) yields for v = v 


j < Z e> 4 | 

(W 

2 

| in j8 

J_ e 

4ttE' 0 

W 

(2kT b ) 3/2 

yp 


-■lA, 



If the particles are, for example, injected at 10 times the average field energy then 


1 r^r 2 - 30 tT 

— mv = — Kl\ 

2 0 2 b 


or 
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II I 


and for deuterons, for example, 



% log/3 


0.773x10" 19 


(Fll) 


This result is used in figure 6. 

It is interesting to compare equation (F9) with equation (33) of reference 8. The 
result of reference 8, derived strictly from a random walk, is for only one absorbing 
wall. The loss rate predicted by equation (F9) (for two absorbing walls) is just twice 
that of reference 7. 

To determine the marginal distribution in 9 for use on figure 7(a) substitute equa- 
tions (F7) and (F8) into (F6). This yields 






when 9 £ 9 s — 

c 2 


~\ 


y (i7) 


when -< 9 < (it - 9 ) 
2 c 


J 


70 



APPENDIX G 


SELECTION OF RANDOM NUMBERS FROM THE NONUNIFORM 
INJECTION DISTRIBUTION OF REFERENCE 5 

The simulation procedure of sampling from a nonuniform distribution when the com- 
puter library contains only a uniformly distributed set of random numbers is quite com- 
mon (refs. 9, 21, and pp. 252-264 of ref. 15). The relation between a uniformly dis- 
tributed random number R and an arbitrarily distributed random number v is 

/ R 

dr = / f(v)dv 

u J 0 

where it is necessary only that the second integral be a monotone increasing function of 
if. This says that the cumulative distribution function R of the probability density f(v) 
is uniformly distributed in Y. If f(v) is integrable, R(^) can be found. But to find 
Y(R) explicitly in the case of interest herein required a root finding method such as, for 
example, the Newton- Raphson method (ref. 22). 

The injection distribution of reference 5 is 


R(Y) = 


C * O -mv 2 /2kT, 9 

I ( (A9) 2 ) e b v 2 



<(A0r>e 


-mv^/2kT 


b v 2 


dv 


where ((A0) 2 ) is given by equation (11c). Letting x= V m /2kT^ v results in a new 
expression for R(y / ): 
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R(*0 = 


/ 


^/m/2kT^ -jT 


r±u 

l X 2 U 


[Xe- x2 *( 

x-i-Wxle- x2<JX 

1_Vj t ' 

v, 2 x/ J 


\ I 2 

- x + (x-l 
V 2x 

> rfX Je' X * 


This was integrated by the method of Gaussian quadratures. The curve fit 


R(f) = erf [1.073 



is a good approximation to the result. This is shown on figure 12 where V = Y/\ 



Figure 12. - Random number distribution to fit injection distribution of 
reference 5. 


1 


used as the abscissa and V m /2kT^ v Q = was set equal to V3/2 to make ref- 

erence velocity v Q equal to the mean field velocity. 

For an initial approximation V(R) for use in the Newton- Raphson method, the curve 

o 

V = arcsin R 


was used. 

Results of this procedure to generate the injection distribution is shown to be satis- 
factory in figure 8(c). 
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